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Abstract 

In  order  to  enhance  the  accuracy  of  the  Global  Positioning  System  (GPS)  it  is 
proposed  to  improve  the  GPS  position  determination  algorithm  used  in  the  GPS  re¬ 
ceiver  to  solve  the  trilateration  equations.  The  use  of  direct,  or  closed-form,  solutions 
based  on  stochastic  modeling  and  estimation  techniques  is  investigated  to  achieve 
this  improvement.  The  objective  of  this  thesis  research  is  the  development  and  exper¬ 
imental  analysis  of  new  closed-form  position  determination  algorithms  that  work  in 
the  presence  of  pseudorange  measurement  noise.  The  mathematical  derivation  of  two 
closed-form  solution  algorithms,  based  entirely  on  linear  mathematics,  is  presented. 
The  closed-form  algorithms  provide  an  estimate  of  the  GPS  solution  parameters 
(viz.,  the  user  position  and  the  user  clock  bias)  as  well  as  the  associated  parameter 
estimation  error  covaxiance.  The  experimental  results  are  based  on  5000  Monte  Carlo 
runs  and  are  produced  through  realistic  simulations  of  typical  NAVSTAR  GPS  low 
noise,  near-earth  navigation  scenarios  and  of  ground-based  pseudolite  planar  array 
scenarios.  Analysis  of  the  experimental  results  is  accomplished  through  direct  com¬ 
parison  to  the  baseline  results  from  the  conventional  Iterative  Least  Squares  (ILS) 
algorithm,  which  is  currently  used  in  GPS.  To  ensure  fair  comparitive  analysis,  all 
the  algorithms  are  tested  under  identical  noise  and  geometry  conditions.  The  closed- 
form  algorithms,  both  of  which  produce  identical  results,  are  extremely  sensitive  to 
noise  in  typical  NAVSTAR  GPS  scenarios,  making  them  unsuitable  for  stand-alone 
use;  however,  they  perform  very  well  at  estimating  horizontal  position  parameters 
in  ground-based  pseudolite  planar  array  scenarios  where  the  ILS  algorithm  breaks 
down  due  to  poor  geometry.  For  the  near-earth  GPS  scenarios,  the  use  of  an  addi¬ 
tional  nonlinear  measurement  in  a  supplementary  algorithm  is  required  to  refine  the 
solution  produced  by  the  closed-form  algorithm.  Thus,  the  derivation  of  two  supple¬ 
mentary  algorithms  is  presented.  The  first  supplementary  algorithm  is  based  on  a 


IX 


maximum  likelihood  approach  and  the  second  uses  a  Kalman-like  update  approach. 
The  maximum  likelihood  algorithm  is  capable  of  producing  solutions  with  accuracy 
equivalent  to  that  of  the  conventional  ILS  algorithm,  both  in  terms  of  mean  esti¬ 
mation  errors  and  experimental  estimation  parameter  standard  deviations;  however, 
the  capability  to  predict  its  estimation  error  covariance  is  lost.  The  performance 
of  the  Kalman  update  algorithm  is  marginally  degraded,  most  notably  in  its  user 
clock  bias  estimation  errors,  but  it  is  capable  of  predicting  its  estimation  error  co- 
variance.  The  capability  to  predict  the  estimation  error  covariance  is  dependent  on 
how  well  the  closed-form  algorithm  estimates  the  pseudorange  measurement  noise 
variance,  which  in  turn  is  a  function  of  satellite  availability.  The  performance  of 
the  closed-form,  supplemented  by  the  Kalman  update,  algorithm  is  comparable  to 
the  ILS  algorithm’s  performance.  The  advantages  introduced  by  the  developed  al¬ 
gorithm  are;  1)  The  capability  to  estimate  its  estimation  error  covariance,  and  2) 
The  potential  for  computational  efficiency  due  to  the  direct,  closed-form,  nature  of 
the  solution. 
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IMPROVED  MATHEMATICAL  MODELING  FOR  GPS  BASED 

NAVIGATION 


1.  Introduction 

1.1  Background 

The  NAVSTAR  Global  Positioning  System  (GPS)  is  an  extremely  accurate 
navigation  system;  nevertheless,  there  are  errors  associated  with  the  GPS  position 
fix,  leaving  room  for  further  improvements.  The  GPS  is  a  space  based  satellite  radio 
navigation  system  that  provides  three  dimensional  (3-D)  user  positioning  by  solving  a 
set  of  nonlinear  pseudorange  trilateration  equations.  A  common  approach  to  solving 
the  nonlinear  equations  is  to  linearize  the  pseudorange  equations  and  calculate  the 
user  position  iteratively,  starting  with  a  user  provided  initial  position  guess  [27].  For 
near-earth  navigation,  the  center  of  the  earth  is  a  good  initial  guess  which  guarantees 
that  the  iterative  algorithm  will  quickly  converge  towards  the  correct  GPS  solution. 
An  area  of  potential  improvement,  that  has  been  the  topic  of  several  papers  in 
recent  years,  is  the  use  of  closed-form  solutions  to  the  nonlinear  GPS  pseudorange 
equations.  The  use  of  closed-form  or  direct  solutions  has  been  motivated  by  the  fact 
that  direct  solutions  shorten  the  computational  cycle  and  that  they  do  not  depend 
on  an  initial  position  guess.  Closed-form  solution  algorithms  have  been  developed  by 
Bancroft  [2],  Krause  [14],  Abel  and  Chaffee  [1],  Chaffee  and  Abel  [6],  and  Hoshen  [11] 
using  different  approaches. 

The  recent  derivations  of  the  closed-form  solutions  considered  only  the  exactly 
determined  case  in  which  four  of  the  available  pseudorange  measurements  are  used  to 
obtain  exact  positioning  solutions.  The  motivation  for  considering  only  four  satellites 
is  questionable  since,  if  a  five  degree  elevation  mask  angle  is  assumed,  it  has  been 
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shown  that  for  near-earth  navigation,  there  are  at  least  five  satellites  in  view  at  all 
times  and  at  least  seven  satellites  in- view  80  %  of  the  time  [28].  An  elevation  mask 
angle  is  used  to  mask  satellites  close  to  the  horizon.  Although  satellite  vehicles  close 
to  the  horizon  are  in- view  by  line-of-sight,  their  transmitted  data  is  prone  to  a  large 
amount  of  atmospheric  interference  caused  by  the  long  propagation  path  through  the 
ionosphere  and  troposphere.  In  airborne  applications,  satellite  visibility  will  tend  to 
improve  with  altitude,  potentially  increasing  the  number  of  satellite  pseudorange 
measurements  available.  In  addition,  some  receivers  have  access  to  signals  from 
the  GLONASS  constellation,  the  former  USSR’s  satellite-based  global  navigation 
system,  potentially  doubling  satellite  availability. 

Previous  works  on  closed-form  solutions  of  the  GPS  pseudorange  equations 
did  not  make  use  of  the  pseudorange  measurements  from  all  in-view  satellites.  This 
research  differs  from  previous  works  in  that  it  will  treat  an  overdetermined  system, 
making  use  of  all  in-view  satellites.  Furthermore,  the  research  specifically  addresses 
the  issue  of  developing  a  reliable  closed-form  solution  that  works  in  the  presence 
of  noise.  With  the  exception  of  [8],  the  pseudorange  equations  have  been  treated 
as  a  deterministic  set  of  equations.  Benefits  can  be  gained  by  recognizing  that 
pseudorange  measurements  are  noise  corrupted;  hence,  the  stochastic  nature  of  the 
measurements  should  be  reflected  in  the  GPS  pseudorange  equations  from  the  onset 
to  develop  a  stochastically  sound  position  solution  estimate.  Furthermore,  it  would 
allow  for  obtaining  an  estimate  of  the  covariance  of  the  position  estimation  error. 

1.2  Research  Motivation 

The  goal  in  the  design  of  a  navigation  system  is  to  obtain  the  best  possible  po¬ 
sitioning  accuracy  by  eliminating  or  at  least  minimizing  the  impact  of  error  sources. 
The  GPS  system  errors  can  be  attributed  to  seven  basic  sources  of  error:  satellite 
clock  errors,  atmospheric  delays,  group  delay,  ephemeris  errors,  receiver  noise  and 
resolution,  multipath  errors,  and  receiver  vehicle  dynamics  [21].  These  sources  of 
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error  are  briefly  addressed  in  the  discussion  of  the  stochastic  modeling  of  the  pseudo¬ 
range  equations  in  Chapter  II.  Although  not  the  only,  or  most  intuitive,  approach 
to  improving  positioning  accuracy,  improvement  of  the  GPS  position  determination 
algorithm  used  within  the  receiver  to  produce  the  position  estimate  from  the  pseudo¬ 
range  measurements  undoubtedly  impacts  the  achievable  accuracy  of  the  solution. 
Improvements  to  the  algorithm  will  complement  any  other  system  improvements  to 
reduce  errors  in  any  of  the  seven  basic  error  sources.  There  are  numerous  potential 
benefits  that  encourage  improvements  to  the  GPS  position  determination  algorithm; 
a  few  that  provided  the  motivation  to  this  thesis  are  discussed  in  this  section. 

1.2.1  Improved  User  Positioning.  GPS  plays  an  extremely  important  role 
in  positioning  applications,  both  military  and  civilian.  Regardless  of  whether  precise 
positioning  service  (PPS)  or  standard  positioning  service  (SPS)  is  being  used,  an 
improved  GPS  position  determination  algorithm  that  provides  a  closed-form  solution 
to  the  GPS  pseudorange  equations  considering  all  in-view  satellites  can  improve  the 
accuracy  of  the  GPS  position  fix.  The  improvement  is  expected  to  be  achieved  as 
a  result  of  computing  an  exact  solution  to  the  nonlinear  equation  as  opposed  to 
introducing  approximations  by  linearizing  the  equations. 

1.2.2  Test  Range  Enhancement.  The  Submeter  Accuracy  Reference  Sys¬ 
tem  (SARS)  navigation  test  reference  system  being  developed  at  the  746*^  Test 
Squadron  at  Holoman  Air  Force  Base,  is  used  primarily  in  the  flight  testing  of  in¬ 
tegrated  navigation  systems  [20].  With  the  improvement  of  GPS,  the  accuracy  of 
the  integrated  aircraft  navigation  systems  is  also  improving.  Since  the  test  refer¬ 
ence  system  must  provide  much  higher ,  accuracy  than  the  system  being  tested,  the 
importance  of  improved  test  range  accuracy  can  never  be  overemphasized.  The  use 
of  an  improved  GPS  position  estimation  algorithm  can  enhance  the  accuracy  of  the 

SARS. 
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1.2.3  Diversified  Applications  for  GPS.  Current  iterative  GPS  algorithms 
are  guaranteed  to  converge  to  a  correct  position  solution  for  near-earth  users,  initial¬ 
izing  the  iterations  at  the  earth’s  center  and  assuming  a  zero  user  clock  bias  [1].  In 
space  applications,  inverted  GPS  applications  such  as  the  SARS  [25],  [20],  and  un¬ 
conventional  applications  that  involve  the  use  of  both  pseudolites  and  satellites,  the 
lack  of  a  sufficiently  good  initial  guess  may  lead  to  convergence  towards  the  wrong 
solution.  Pseudolite  are  ground-based  transmitters  that  provide  GPS-like  position¬ 
ing  data  and  can  be  used  for  augmentation.  Inverted  GPS  applications,  and  the 
use  of  pseudolites  are  applications  that  are  currently  being  considered  by  the  746*^ 
Test  Squadron  for  the  SARS  as  a  result  of  work  performed  at  AFIT  by  McKay  [20]. 
The  fact  that  a  closed-form  solution  will  not  require  an  initial  position  guess  is  an 
advantage  for  these  applications. 

1.2.4  Position  Estimate  Error  Covariance.  By  stochastically  modeling 
the  GPS  pseudorange  equations,  solving  for  position  becomes  a  stochastic  estima¬ 
tion  problem.  The  use  of  correct  stochastic  modeling  and  of  Kalman  Filtering  like 
techniques  to  solve  the  estimation  problem  will  lead  to  a  GPS  solution  that  pro¬ 
vides  accurate  estimates  of  the  position  estimation  error  covariance  in  addition  to 
the  position  estimate  itself.  This  will  introduce  a  new  confidence  factor  into  GPS 
positioning  and  is  critical  for  the  integration  of  GPS  with  additional  sensors  for  inte¬ 
grated  navigation  systems.  This  will  specifically  enhance  the  accuracy  of  navigation 
systems  that  are  integrated  using  a  federated  approach  but  will  not  help  in  deep  in¬ 
tegration  schemes  as  discussed  in  [16],  where  the  raw  pseudorange  signals  are  being 
used  for  system  integration.  Federated  system  integration  is  done  at  processed  data 
levels  using  computed  positioning  data  from  the  navigation  sensors.  This  approach 
is  common  when  system  integration  is  performed  as  an  upgrade  on  existing  sensors 
that  do  not  provide  access  to  raw  data  signals  required  for  deep  integration  [9]. 
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1.2.5  Computational  Efficiency.  The  use  of  direct,  or  closed-form  solu¬ 
tion  algorithms  tend  to  be  more  computationally  efficient  than  iterative  or  recursive 
algorithms.  This  efficiency  results  in  reduction  of  the  computational  cycle  which 
is  most  advantageous  for  fast  moving  vehicles.  The  use  of  a  positioning  algorithm 
with  short  computation  cycles  is  of  extreme  importance  to  space  vehicles  where  the 
earth’s  rotation  and  the  vehicle  movement,  in  the  computation  interval,  needs  to  be 
taken  into  account. 

1.3  Problem  Statement 

The  objective  of  this  thesis  is  to  provide  a  closed-form  solution  to  the  GPS 
pseudorange  equations  and  present  a  way  of  using  it  effectively  in  the  presence  of 
GPS  measurement  noise.  This  thesis  work  develops  an  improved  closed-form  mathe¬ 
matical  solution  to  the  GPS  pseudorange  equations,  implements  an  algorithm  based 
on  the  mathematical  solution,  and  performs  an  experimental  analysis  of  the  algo¬ 
rithm  using  realistic  Monte  Carlo  simulations.  The  mathematical  derivation  of  the 
closed-form  solution  used  in  this  research  is  closely  based  on  notes  provided  by 
Dr.  Meir  Pachter  [22]. 

The  closed-form  positioning  solution  will  consider  pseudorange  measurements 
from  all  available  satellites,  to  obtain  a  position  fix,  since  they  are  readily  handled 
by  the  closed-form  solution  algorithm.  A  minimum  of  five  satellites  will  be  needed  to 
obtain  an  initial  3-D  position  fix  using  the  algorithm  developed  in  this  thesis,  but  this 
is  not  a  serious  problem  since  there  are  always  more  than  four  satellites  in- view  when 
the  NAVSTAR  GPS  constellation  is  fully  operational.  The  thesis  work  will  emphasize 
the  proper  use  of  stochastic  modeling  and  estimation  in  order  to  provide  appropriate 
weighting  of  satellite  pseudorange  data  in  producing  a  positioning  solution. 
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1-4  Scope 

The  focus  of  this  research  is  on  the  improvement  of  the  GPS  algorithm  to  obtain 
better  positioning  accuracy  and  to  obtain  shorter  computation  cycles.  There  are 
numerous  other  factors  in  GPS  receiver  design  that  affect  GPS  positioning  accuracy 
that  will  not  be  addressed  in  this  thesis.  The  GPS  will  be  considered  at  the  system 
level,  hence  the  inner  workings  of  the  GPS  receiver  will  not  be  considered.  The 
GPS  pseudoranges  that  will  be  used  in  the  research  are  corrected  pseudoranges.  The 
corrected  pseudoranges  will  represent  the  pseudoranges  as  they  would  be  provided  by 
the  receiver  after  all  known  correction  factors  have  been  applied  and  known  errors 
modeled  out  of  the  raw  pseudorange  measurements.  The  corrections  of  these  errors 
are  GPS  receiver  design  issues  beyond  the  scope  of  this  thesis.  Throughout  this 
thesis,  the  term  pseudoranges  is  treated  as  meaning  corrected  pseudoranges. 

The  experimental  analysis  will  be  limited  to  obtaining  positioning  performance 
indexes  relative  to  the  conventional  ILS  GPS  algorithm.  The  experimental  analysis 
will  not  encompass  any  computational  and  numerical  performance  issues  that  would 
have  to  be  addressed  prior  to  using  the  algorithm  in  real-time  applications.  Due 
to  the  proprietary  nature  of  the  algorithms  in  current  commercially  available  GPS 
receivers,  comparisons  to  actual  GPS  receiver  performance  will  not  be  attempted. 

1.5  Assumptions 

Assumptions  must  be  made  about  the  noise  corrupting  the  pseudorange  mea¬ 
surements  in  order  to  allow  the  use  of  a  simple  stochastic  model  and  to  simplify  the 
stochastic  estimation  problem.  The  following  assumptions  are  used  in  this  thesis: 

1.  After  all  known  corrections  are  applied  to  the  pseudorange  measurements,  the 
residual  noise  corrupting  the  pseudorange  measurement  is  a  zero-mean  Gaus¬ 
sian  distributed  noise. 

2.  The  noises  on  all  the  pseudorange  measurements  have  equal  variance. 
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3.  The  noises  on  the  pseudorange  measurements  are  uncorrelated  with  one  an¬ 
other. 

4.  The  effects  of  Selective  Availability  (SA)  on  pseudorange  noise  is  not  consid¬ 
ered. 

5.  A  10  degree  elevation  angle  is  always  achievable  for  determination  of  in-view 
GPS  satellites. 

6.  The  GPS  satellite  constellation  has  24  fully  operational  satellites. 

The  extent  to  which  the  assumptions  are  valid  is  not  exactly  known;  con¬ 
sequently,  the  impact  of  using  these  assumptions  can’t  be  determined  beforehand. 
Some  of  the  assumptions  made  are  necessary  to  obtain  the  solution  to  the  pseudo¬ 
range  equations.  Others  are  required  to  establish  a  realistic  baseline  for  satellite 
availability  used  for  the  experimental  simulations.  Attempts  will  be  made  to  qualify 
the  significance  of  these  assumptions  on  the  positioning  solution  through  experimen¬ 
tal  analysis. 

1 . 6  Methodology 

This  thesis  contains  three  distinct  phases.  The  first  phase  is  the  development 
of  the  closed-form  solution  to  the  pseudorange  equations.  The  second  phase  is  the 
implementation  of  the  experimental  algorithms  based  on  the  developed  solutions. 
The  third  phase  is  the  experimental  analysis  which  can  only  be  initiated  after  the 
algorithms  in  the  second  phase  are  fully  completed,  implemented,  and  debugged. 
The  approach  will  be  iterative  in  nature  since  rework  of  the  mathematical  solution 
may  be  required  after  some  experimental  analysis  is  performed,  which  in  turn  will 
require  changes  to  the  algorithms  and  the  performing  of  new  simulations. 

1.6.1  Phase  1  -  Development  of  a  Closed  Form  Solution.  This  phase  con¬ 
sists  of  the  mathematical  derivation  of  a  closed-form  solution  to  the  GPS  pseudorange 
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equations.  The  solution  will  consider  n  >  5  satellites  and  will  capture  the  stochastic 
nature  of  the  pseudorange  measurement. 

The  corrected  pseudorange  can  be  modeled  as  the  true  Euclidean  range  with  an 
unknown  clock  bias  and  some  measurement  noise  superimposed;  thus,  the  stochastic 
pseudorange  measurement  equation  is  given  by 


Ri  =  yl{ux-  XiY  +  {uy  -  yif  +  (uz  -  ZiY  +  b  +  Wi,  (1.1) 

where  Ri  is  the  corrected  pseudorange  corresponding  to  the  satellite,  {ua;,Uy,Uz) 
are  the  user  position  coordinates,  {xi,  yi,  Zi)  are  the  known  coordinates  of  the  satel¬ 
lite,  b  is  the  range-equivalent  user  clock  bias  defined  as  the  error  in  the  receiver  clock 
times  the  speed  of  light,  and  Wi  is  zero-mean  Gaussian  pseudorange  measurement 
noise.  It  is  reasonable  to  assume  that  all  measurements  are  subject  to  the  same 
noise;  therefore,  they  will  have  the  same  variance,  <7^.  However,  the  noise  terms 
are  not  correlated  between  satellites.  The  point  of  departure  for  the  mathematical 
derivation  of  the  closed-form  stochastic  solution  is  Equation  (1.1). 

One  such  measurement  equation  is  available  for  each  of  the  n  in-view  satel¬ 
lites.  The  n  equations  will  be  referred  to  as  the  GPS  pseudorange  equations.  All 
positions  used  in  the  derivation  will  be  expressed  in  Earth  Centered  Earth  Fixed 
(ECEF)  coordinates.  These  n  nonlinear  equations  will  be  solved  algebraically  for 
the  estimated  user  position,  (ux,Uy,Uz)  and  user  clock  bias,  b.  This  will  be  achieved 
through  algebraic  manipulation  to  reduce  the  GPS  pseudorange  equations  into  a 
linear  regression  in  the  form  of  the  standard  linear  measurement  model  as  defined 
in  [18].  The  linear  regression  is  given  by 

Z  =  HX  +  V,  (1.2) 

where  Z  is  the  measurement  vector,  X  is  the  vector  of  unknown  quantities  com¬ 
prising  the  user  position  coordinates  (ux,Uy,Uz)  and  the  user  clock  bias  (6),  H  is 
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the  regressor,  an  (n  x  4)  measurement  matrix,  and  V  is  a  Gaussian  noise  vector,  or 
equation  error  vector,  whose  covariance  matrix  must  be  determined.  The  solution  is 
obtained  by  solving  the  static  estimation  problem.  The  covariance  of  the  estimation 
error  of  the  user  position  and  user  clock  bias  is  also  obtained. 

1.6.2  Phase  2  -  Implementation  of  the  Algorithms.  Matlab  [17]  algorithms 
based  on  the  solution  processes  to  the  GPS  pseudorange  equations  developed  in 
Phase  1  are  implemented.  The  main  algorithm  will  be  referred  to  as  the  closed-form 
algorithm  to  distinguish  it  from  conventional  iterative  algorithms.  The  closed-form 
algorithm  will  expect  as  inputs  pseudorange  measurements  from  n  >  5  satellites 
and  will  provide  as  a  solution  the  estimates  of  the  user  position  coordinates,  user 
clock  bias,  and  the  associated  estimation  error  covariance  matrix.  The  closed-form 
algorithm,  as  implied  by  its  name,  will  be  non-iterative  in  nature;  hence,  will  not  be 
using  an  initial  guess  of  the  receiver  position  and  clock  bias. 

Two  additional  algorithms  based  on  conventional  iterative  solutions  to  the  GPS 
pseudorange  equations  will  also  be  implemented  in  Matlab.  Both  algorithms  will  be 
iterative  in  nature  and  will  both  use  the  linearized  pseudorange  equation  approach 
described  in  [27].  These  algorithms  will  expect  the  same  input  that  was  provided 
to  the  direct  algorithm,  the  pseudorange  measurements  from  n  >  5  satellites.  Since 
these  algorithms  require  an  initial  guess  of  the  receiver  position  and  clock  bias,  the 
initial  guess  will  be  built  into  the  algorithms.  Due  to  the  lack  of  a  better  guess  it  is 
common  practice  to  use  the  center  of  the  earth  as  the  initial  position  guess  and  zero 
as  the  initial  clock  bias  guess. 

The  first  iterative  algorithm  will  use  all  available  pseudorange  measurements 
to  calculate  a  least  squares  solution;  the  algorithm  will  be  referred  to  as  the  Iterative 
Least  Squares  (ILS)  algorithm.  The  second  algorithm  will  provide  a  solution  based  on 
only  the  best  four  pseudorange  measurements  based  on  geometry  strength.  Dilution 
of  Precision  (DOP)  is  used  as  a  measure  of  geometrical  strength  [5];  thus,  a  frequently 
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used  approach  is  to  select  the  four  satellites  whose  geometry  yields  the  smallest 
Position  Dilution  of  Precision  (PDOP),  which  considers  only  the  effect  of  the  three 
ECEF  axis  position  error, s  [27].  To  simplify  the  experimental  data  collection,  the 

■I.. 

Geometric  Dilution  of  Precision  (GDOP)  will  be  used  to  evaluate  the  geometry 
since  GDOP  will  be  the  same  as  PDOP  if  a  zero  receiver  clock  drift  rate  is  assumed. 
GDOP  is  defined  as, 

GDOP  =  (1.3) 

where  H  is  the  matrix  of  direction  cosines  associated  with  the  conventional  iterative 
algorithms.  The  algorithm  which  considers  only  the  four  satellites  which  yield  the 
lowest  GDOP  is  referred  to  as  the  Iterative  Best  Four  (IBF)  algorithm.  Both  algo¬ 
rithms  provide  a  solution  that  includes  the  estimated  receiver  position  coordinates, 
and  the  receiver  clock  bias.  These  algorithms  do  not  have  the  capability  of  providing 
an  indication  of  how  good  the  calculated  estimates  actually  are  unless  they  are  arti¬ 
ficially  informed  of  the  pseudorange  noise  strength.  Since  in  a  simulated  experiment 
the  noise  strength  is  known,  it  is  possible  to  inform  the  algorithm  of  the  value  of 
(7  to  verify  how  well  the  algorithm’s  predicted  performance  compares  to  its  actual 
performance  based  on  the  experimentally  obtained  covariance  yielded  over  the  5000 
Monte  Carlo  runs.  A  discussion  of  how  the  number  of  Monte  Carlo  runs  was  selected 
is  presented  in  Chapter  IV. 

( 

1.6.3  Phase  3  -  Experimental  Analysis.  The  experimental  analysis  phase 
will  contain  the  bulk  of  the  work  of  this  study.  The  object  of  this  phase  is  to 
provide  thorough  and  realistic  performance  analyses  of  the  closed-form  algorithm 
implemented  in  the  second  phase.  In  addition  to  evaluating  how  good  the  position 
and  clock  bias  estimates  are,  particular  emphasis  will  be  placed  on  how  well  the  algo¬ 
rithm  believes  it  is  doing.  The  provision  of  the  estimation  error  covariance  is  the  key 
innovation  to  the  algorithm  developed  in  this  thesis;  hence,  it  is  of  extreme  impor¬ 
tance  that  the  validity  of  this  feature  be  established.  Some  performance  comparisons 
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between  the  direct  algorithm  and  the  two  iterative  algorithms  will  be  performed  to 
provide  an  appreciation  of  the  relative  performance. 

GPS  ephemeris  data  will  be  generated  using  GPS  simulation  tools  available  in 
the  GPSoft  Satellite  Navigation  Toolbox  for  Matlab  [10].  This  will  provide  realistic 
data  from  which  true  ranges  from  all  in-view  GPS  satellites  to  an  arbitrarily  selected 
receiver  position  can  be  calculated.  After  adding  an  arbitrary  clock  bias  to  all  the 
ranges,  a  zero-mean  random  noise  of  preselected  variance  is  superimposed  to  repre¬ 
sent  the  Gaussian  measurement  noise.  These  generated  pseudoranges  will  serve  to 
perform  the  analysis  of  the  GPS  algorithm  through  the  running  of  5000  Monte  Carlo 
trials.  Using  the  Satellite  Navigation  Toolbox,  the  capability  exists  of  simulating 
realistic  noise  corrupted  pseudorange  measurements  which  could  be  applied  directly 
to  the  GPS  position  determination  algorithm  as  would  be  the  case  in  a  real  world 
scenario.  The  approach  used  of  simulating  just  the  GPS  satellite  ephemeris  data  and 
producing  the  simulated  pseudoranges  was  preferable  for  the  following  reasons: 

1.  It  provides  a  more  structured  data  set  for  analysis  of  the  algorithms  since  only 
the  desired  effects  are  being  considered  and  the  amount  of  corruption  on  the 
pseudorange  measurements  is  exactly  known;  and 

2.  Since  the  pseudoranges  are  produced  starting  from  exactly  known  position 
coordinates,  comparisons  against  the  true  position  for  determining  accuracy 
axe  possible. 

Numerous  trials  will  be  run  under  different  scenarios  to  provide  a  thorough 
picture  of  the  performance.  As  an  end  result  of  the  Monte  Carlo  simulations  in 
Matlab,  an  estimate  of  the  position  fix  accuracy  attainable  using  the  developed 
algorithm  will  be  obtained  and  consequently  compared  to  that  achievable  using  the 
conventional  algorithms.  Recognizing  that  a  large  number  of  Monte  Carlo  trials  are 
required  to  achieve  statistically  reliable  experimental  results,  all  simulations  will  be 
based  on  5000  Monte  Carlo  runs. 
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Through  the  experimental  trials  the  issue  of  geometry  will  also  be  analyzed. 
Experimental  scenarios  with  varying  degrees  of  GDOP  will  be  used  for  the  trials 
to  bring  out  the  effects  of  geometry  on  GPS  positioning.  GDOP  is  commonly  used 
as  a  measure  of  GPS  satellite  geometry,  whereby  a  lower  GDOP  is  an  indication 
of  better  geometry;  hence,  lower  positioning  error.  As  suggested  by  McKay,  GDOP 
may  not  be  the  best  measure  of  geometry  [20] ;  therefore,  the  condition  number  of  the 
regressor  matrix  will  also  be  examined  and  attempts  will  be  made  to  bring  forth  any 
relationships  to  the  accuracy  bounds.  It  is  expected  that,  through  the  investigation 
of  satellite  geometry,  sorne  interesting  issues  related  to  the  accuracy  of  the  GPS 
positioning  calculations  will  arise.  The  validity  of  considering  only  the  four  satellites 
that  yield  the  best  GDOP  to  calculate  a  position  fix,  as  well  as  the  validity  of  using 
GDOP  as  the  measure  of  geometry,  will  be  examined. 

1.7  Thesis  Overview 

Chapter  II  presents  the  background  theory  upon  which  this  research  will  be 
based.  The  first  part  includes  a  thorough  literature  review  that  provides  a  summary 
of  current  knowledge  in  the  field  of  closed-form  GPS  algorithms  and  the  application 
of  stochastic  modeling  to  the  GPS  pseudorange  equations.  The  second  part  discusses 
the  basic  concepts  of  GPS  technology.  Particular  emphasis  is  placed  on  the  formula¬ 
tion  of  the  GPS  pseudorange  equations  and  the  theory  behind  the  current  iterative 
algorithms. 

Chapter  III  presents  a  thorough  development  of  the  closed-form  mathematical 
solution  to  the  GPS  pseudorange  equations.  The  development  of  supplementary  al¬ 
gorithms  that  make  use  of  a  nonlinear  measurement  equation  and  that  were  deemed 
essential  to  enhance  the  performance  of  the  closed-form  algorithm  is  also  addressed. 
The  purpose  of  these  enhanced  algorithms  is  to  produce  an  improvement  to  the  solu¬ 
tion  achieved  by  the  closed-form  algorithm;  they  are  not  capable  of  providing,  nor  are 
they  intended  to  provide,  stand-alone  solutions.  The  two  stage  algorithms  that  are 
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presented  are  the  Maximum  Likelihood  approach  and  the  Kalman  Update  approach. 
The  background  theory  and  the  operating  concept  behind  each  approach  is  discussed 
and  the  complete  mathematical  derivations  are  presented.  The  development  of  an 
alternate  closed-form  algorithm  is  also  presented. 

Chapter  IV  presents  the  experimental  portion  of  the  thesis.  The  first  part 
discusses  how  the  experiment  was  set  up  and  how  the  Monte  Carlo  trials  were  run. 
Next  is  presented  a  brief  discussion  of  how  the  results  are  intended  to  be  interpreted. 
The  experimental  results  are  then  presented  and  the  chapter  sums  up  with  a  detailed 
analysis  of  the  results. 

Chapter  V  presents  a  brief  summary  of  the  performance  related  issues  of  the 
closed-form  algorithm  and  the  supporting  algorithms  developed  for  this  thesis  re¬ 
search.  Emphasis  is  placed  on  the  identifying  the  areas  of  strength  for  the  algorithms 
and  suggesting  applications  for  which  they  are  best  suited.  The  shortcomings  seen 
in  the  Chapter  IV  related  to  algorithm  deficiency  and  inconclusive  findings  are  dis¬ 
cussed,  and  potential  approaches  to  correcting  the  problems  are  presented.  The 
chapter  sums  up  with  recomendations  for 
bringing  the  concept  of  using  a  stochastic 
algorithm  closer  to  reality. 


future  work  that  can  prove  rewarding  in 
closed-form  GPS  position  determination 
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11.  Background 

This  chapter  presents  the  background  theory  upon  which  this  research  will  be 
based.  The  first  part  includes  a  thorough  literature  review  that  provides  a  summary 
of  current  knowledge  in  the  field  of  closed-form  GPS  algorithms  and  the  application 
of  stochastic  modeling  to  the  GPS  pseudorange  equations.  The  second  part  discusses 
the  basic  concepts  of  GPS  technology.  Particular  emphasis  is  placed  on  the  formula¬ 
tion  of  the  GPS  pseudorange  equations  and  the  theory  behind  the  current  iterative 
algorithms. 

2.1  The  Global  Positioning  System 

This  section  presents  an  overview  of  the  Global  Positioning  System  to  provide 
some  insight  into  the  complexity  of  the  system.  This  section  serves  to  focus  atten¬ 
tion  to  the  specific  portions  of  GPS  of  interest  in  this  thesis,  namely  the  position 
determination  algorithm  within  the  GPS  receiver,  and  shows  how  it  fits  into  context 
of  the  overall  GPS  system.  GPS  specific  terminology  and  the  coordinate  systems 
used  in  this  thesis  are  discussed  as  well. 

2.1.1  GPS  System  Overview.  GPS  is  a  satellite  based  radio-navigation 
system  that  provides  worldwide,  virtually  continuous,  three-dimensional  (3-D)  posi¬ 
tioning  and  accurate  timing.  The  beauty  of  the  system  is  in  its  apparent  simplicity 
since,  from  the  user’s  perspective,  extremely  accurate  positioning  can  be  achieved 
with  the  use  of  a  simple,  fairly  inexpensive  GPS  receiver.  For  these  reasons,  GPS 
is  rapidly  becoming  the  positioning  sensor  of  choice  for  both  military  and  civilian 
users.  The  continuous  worldwide  coverage  provided  by  GPS  makes  it  ideally  suited 
for  air,  land,  and  sea  navigation  applications.  It  must  be  recognized  that  there  is 
much  more  to  GPS  than  the  portion  that  the  typical  user  is  dealing  with  in  obtain¬ 
ing  a  GPS  position  fix.  GPS  is  composed  of  three  segments:  the  space  segment, 
the  control  segment,  and  the  user  segment.  Each  segment  is  essential  to  the  proper 
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functioning  of  GPS  as  an  accurate  and  reliable  navigation  tool.  The  typical  GPS 
user  is  concerned  with  the  user  segment  only. 

The  GPS  space  segment  refers  to  the  GPS  satellite  constellation.  The  GPS 
constellation  consists  of  24  satellites,  21  active  satellites  plus  three  active  spares,  in 
six  orbital  planes.  Each  of  the  satellites  has  an  orbital  period  of  approximately  12 
hours  (half  a  sidereal  day).  A  healthy  GPS  constellation  provides  satellite  coverage 
such  that,  for  near-earth  locations,  there  are  always  at  least  five  satellites  in- view  and 
at  least  seven  satellites  in- view  80  percent  of  the  time  [28].  The  satellites  transmit 
time-tagged  navigation  messages  which  the  GPS  receivers  (user  segment)  use  to  cal¬ 
culate  their  positions.  The  navigation  message  information  required  by  the  receivers 
to  perform  their  function,  includes  GPS  time,  satellite  ephemeris  data,  correction 
data,  and  system  almanac  data  [27]. 

The  control  segment  is  composed  of  five  stations  spread  over  the  world,  which 
monitor  and  control  the  satellite  orbits  and  GPS  time.  Only  one  of  the  five  sta¬ 
tions  is  the  master  control  station  and  only  three  of  the  remaining  four  stations  are 
uplink  stations  capable  of  transmitting  data  back  to  the  satellites.  The  five  sta¬ 
tions  receive  the  same  signals  seen  by  all  users  and  collect  pseudoranges  to  all  the 
satellites.  All  pseudorange  measurements  collected  by  the  stations  are  transmitted 
to  the  master  control  station  which  then  computes  the  true  satellite  positions  and 
true  GPS  time.  This  is  possible  since  the  stations  are  situated  at  very  well  surveyed 
positions.  The  master  control  station  then  calculates  corrections  for  the  satellites 
and  transmits  them  to  the  uplink  stations  where  the  corrections  are  transmitted  to 
the  satellites  [26]. 

The  user  segment  is  the  GPS  receiver.  The  receiver  receives  the  navigation 
message,  extracts  the  data  and  applies  the  corrections  to  obtain  the  pseudorange 
and  pseudorange  rate  measurements.  The  user  position,  velocity,  and  time  are  then 
obtained  through  an  algorithm  that  calculates  a  solution  from  the  corrected  pseudo- 
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range  measurements.  The  proposed  research  will  specifically  address  an  algorithm 
that  can  be  used  by  a  GPS  receiver  to  obtain  a  position  fix. 

2.1.2  The  Earth  Centered  Earth  Fixed  (ECEF)  Coordinate  Frame.  The 
ECEF  coordinate  frame  is  an  orthogonal  frame  with  its  origin  at  the  earth’s  center 
of  mass.  The  ECEF  frame  is  fixed  to  the  earth  ajid  therefore  rotates  with  the  earth. 
This  frame  consists  of  three  axes:  x,  y,  and  z.  The  z  axis  is  aligned  with  the  earth’s 
spin  axis  directed  north  and  the  x  and  y  axes  lie  in  the  equatorial  plane.  The  x  axis  is 
directed  through  the  Greenwich  Meridian  (0  longitude)  and  the  y  axis  through  the  90 
east  longitude  [27].  ECEF  coordinatization  is  commonly  used  in  GPS  since,  in  near- 
earth  navigation,  the  navigator  wants  to  know  his  positioning  with  respect  to  the 
earth.  Calculations  in  the  GPS  receiver  are  normally  performed  in  the  ECEF  frame 
for  convenience  but  are  converted,  in  the  GPS  receiver,  to  a  coordinate  system  the 
user  selects  for  display.  Geographic  coordinates  (Latitude,  Longitude  and  Altitude) 
are  commonly  used  for  display  purposes. 

2.2  Summary  of  Current  Knowledge 

The  proposed  thesis  research  involves  the  development  and  evaluation  of  a 
closed-form  solution  to  the  GPS  pseudorange  equations  using  stochastic  modeling 
and  estimation.  A  thorough  literature  review  was  required  to  identify  areas  and 
approaches  that  have  not  yet  been  explored  in  the  field  of  closed-form  GPS  solutions 
and  to  establish  the  framework  for  this  thesis  research.  This  review  of  current  liter¬ 
ature  on  the  NAVSTAR  GPS  emphasizes  the  development  of  closed-form  solutions 
to  the  pseudorange  equations  and  stochastic  estimation. 

Relevant  papers  related  to  GPS  closed-form  solutions  published  over  the  past 
twelve  years,  past  theses,  articles,  and  books  covering  more  general  aspects  of  GPS, 
were  reviewed.  The  information  collected  is  critical  to  establishing  that  a  problem 
exists  and  that  the  proposed  thesis  research  is  a  potential  solution  for  improved 
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mathematical  modeling  for  GPS.  The  literature  review  covers  the  following  areas: 
GPS  Overview,  GPS  equations,  conventional  GPS  positioning  solutions,  and  recent 
alternate  approaches  to  GPS  positioning  solutions.  Upon  completion  of  the  literature 
review,  areas  of  interest  that  have  not  yet  been  addressed  will  be  identified  and  the 
viability  of  pursuing  the  thesis  research  will  be  confirmed. 

2.SJ  The  GPS  Pseudorange  Equations.  GPS  uses  radio  signal  timing  to 
measure  the  ranges  between  the  satellites  and  the  GPS  receiver,  making  it  a  time- 
difference-of-arrival  system.  If  ranges  were  being  measured  directly,  we  would  be 
dealing  with  a  multilateration  system  described  by  a  sphericcil  geometry  and  obtain¬ 
ing  a  position  fix  would  be  trivial  [6].  Under  ideal  error  and  noise-free  conditions,  if 
both  the  satellite  and  the  GPS  receiver’s  clock  were  perfectly  synchronized  on  GPS 
time  with  no  error,  then  the  measured  range  would  be  the  true  range  [29].  The  GPS 
receiver  actually  measures  pseudoranges  which  are  corrupted  by  the  receiver  clock 
bias,  noise,  and  other  error  sources.  The  errors  include  atmospheric  delays,  satellite 
clock  errors,  ephemeris  errors,  and  receiver  induced  errors.  The  receiver  clock  bias 
caused  by  the  difference  between  the  receiver  clock  time  and  GPS  time  is  by  far 
the  largest  contributor  to  the  difference  between  pseudorange  and  range.  However, 
the  receiver  clock  bias  is  common  to  a  set  of  simultaneous  measurements,  making  it 
possible  to  treat  it  as  an  unknown  variable  to  be  estimated  along  with  the  receiver 
position,  hence  the  GPS  solution  considers  space  and  time. 

According  to  Parkinson  [24],  the  GPS  pseudorange  equation  that  reflects  all 
the  known  sources  of  error  is  given  by 

Ri  =  {u^  -  XiY  -b  {uy  -  yiY  -I-  -  ZiY  +  c{hu  -  Bi)  -f  c{T  1)  +  E  -\-W,  (2.1) 

where  Ri  is  the  raw  pseudorange  from  the  user  to  the  satellite,  {ux,Uy,Uz)  are 
the  user  position  ECEF  coordinates,  {xi,yi,Zi)  are  the  ECEF  coordinates  of  the 
satellite  as  calculated  from  the  Keplerian  parameters  in  the  satellite’s  ephemeris  data. 
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c  is  the  speed  of  light  in  vacuum,  bu  is  the  receiver  clock  bias,  Bi  is  the  error  in  satellite 
time,  T  is  the  tropospheric  delay,  I  is  ionospheric  delay,  E  are  ephemeris  errors,  and 
W  represents  other  errors  that  can  attributed  to  the  receiver  including  receiver  noise, 
code  loop  quantization  error,  multipath  effects,  and  interchannel  errors. 

Ephemeris  corrections  provided  to  the  satellites  from  the  control  segment  can 
be  used  to  eliminate  the  satellite  time  error  and  the  ephemeris  errors  partially. 
Known  tropospheric  and  ionospheric  error  model  corrections  can  be  applied  to  com¬ 
pensate  for  tropospheric  and  ionospheric  delay  errors  partially.  Improved  receiver 
design  techniques  are  used  to  minimize  the  effects  of  the  receiver  related  errors, 
including  multipath  errors.  Given  the  current  state  of  GPS  receiver  design  technol¬ 
ogy,  the  residual  errors  that  remain  uncompensated  can  be  assumed  to  be  negligibly 
small.  Furthermore,  if  the  residual  errors  are  grouped  together  under  one  random 
variable  w  the  equation  reduces  to  the  expression  presented  as  Equation  (1.1)  known 
as  the  GPS  pseudorange  equation.  If  the  residual  errors  are  neglected  entirely,  the 
ideal  GPS  pseudorange  equation  can  be  expressed  as: 

Ri  «  {u^  -  XiY  -t-  {Uy  -  yiY  -1-  {Uz  -  ZiY  -I-  b.  (2.2) 

The  ideal  GPS  pseudorange  equation  is  a  nonlinear  equation  in  four  unknowns, 
the  three  receiver  position  coordinates,  (uj,,  Uy^u^)  and  the  receiver  clock  bias,  6.  This 
equation  is  the  basis  for  deriving  the  conventional  iterative  GPS  position  solutions. 
At  least  four  GPS  pseudorange  equations  are  required  to  calculate  the  four  unknowns, 
(ux,Uy,U;i,b)  and  obtain  a  solution. 

2.2.2  Conventional  Rerative  Solution  to  the  GPS  Pseudorange  Equations. 
The  conventional  approach  to  solving  the  GPS  pseudorange  equations  is  to  linearize 
Equation  (2.2)  about  a  nominal  solution  for  the  vector  of  unknowns.  The  vector  of 
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unknowns,  and  its  associated  nominal  values,  are  defined  as: 


X  =  [Uj;  Uy  Uz  bf 


and 


“^zn  bn] 


respectively. 

As  described  in  [27],  performing  a  Taylor  series  expansion  of  the  GPS  pseudo¬ 
range  equation  and  ignoring  the  second  and  higher  order  terms,  the  following  equa¬ 
tion  is  obtained: 


Ri  ^  Rjii  “h  ~  Aujc  +  —  —/^Uy  +  —  H"  A6,  (2.3) 
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Considering  the  exactly  determined  case  of  four  pseudorange  equations,  the 
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The  Uij  entries  in  the  H  matrix  are  recognized  as  the  cosine  of  the  angle  between 
the  line-of-sight  vector  from  the  receiver  to  the  satellite  and  the  axis  of  the 
ECEF  coordinate  frame  [27]. 

The  linearized  pseudorange  equations  can  be  written  in  a  more  compact  form 
as  HAX  =  AR.  Solving  for  AX  gives  the  result  AX  =  H~^AR.  This  equation  can 
be  solved  iteratively  using  the  following  procedure: 

1.  Estimate  an  initial  a  nominal  receiver  position  and  clock  bias. 

2.  Calculate  nominal  pseudoranges  and  difference  them  with  measured  pseudo¬ 
ranges  to  obtain  AR. 

3.  Compute  direction  cosines  to  form  the  H  matrix. 

4.  Compute  AX  =  H~^AR. 

5.  Add  AX  to  Xn  forming  a  new  corrected  X„  and  go  back  to  step  2. 

6.  Continue  process  until  convergence  to  a  solution  is  achieved  by  verifying  that 
II  AX II  fti  0  or  that  an  established  threshold  is  attained. 

Upon  completion,  X„,  the  nominal  position  and  bias,  represents  the  best  esti¬ 
mate  of  receiver  position  and  receiver  clock  bias.  This  method  converges  to  a  solution 
within  three  to  five  iterations  even  when  the  initial  position  guess  is  nowhere  close 
to  the  true  position  [14],  for  instance  the  center  of  the  earth.  Drawbacks  of  using 
this  iterative  approach  include  the  approximative  nature  of  the  linearized  equations, 
computational  loading  associated  with  the  inversion  of  a  four  by  four  matrix,  the  re¬ 
quirement  for  an  initial  guess,  and  the  possibility  of  converging  to  the  wrong  solution 
if  the  initial  position  guess  were  not  sufficiently  close  to  the  true  position  [1]  [14]. 
The  last  concern  is  not  an  issue  for  neax-earth  navigation  since  a  unique  solution  is 
guaranteed  if  the  earth’s  center  is  used  as  the  initial  position  guess  and  a  zero  initial 
clock  bias  is  assumed,  but  is  a  serious  concern  if  the  receiver  position  is  outside  the 
GPS  satellite  constellation  where  a  unique  solution  is  not  guaranteed.  The  same 


2-7 


applies  to  certain  inverted  GPS  arrangements.  To  alleviate  such  concerns,  direct 
closed-form  solutions  to  the  GPS  pseudorange  equations  are  sought. 

S.2.3  Closed-Form  Solutions  to  the  GPS  Pseudorange  Equations.  Although 
closed-form  solutions  to  the  GPS  pseudorange  equations  are  attractive,  the  concept  is 
not  new.  Joseph  Hoshen  [11]  proposed  that  a  closed-form  solution  to  two-dimensional 
equations  in  the  form  of  the  GPS  pseudorange  equations  may  have  been  available 
since  the  third  century  BC  in  the  form  of  the  Problem  of  Apollonius.  Since  GPS  is 
a  fairly  recent  system,  the  first  article  in  the  open  literature  concerned  with  closed- 
form  solutions  specifically  tailored  to  the  GPS  pseudorange  equations  is  Stephen 
Bancroft’s  in  1985  [2].  Bancroft  [2]  developed  an  algebraic  solution  to  the  GPS 
pseudorange  equations  that  was  noniterative  in  nature.  His  method  provides  an  ex¬ 
act  solution  in  the  exactly  determined  system  using  four  satellites;  like  the  iterative 
solution,  however,  it  provides  a  least  squares  solution  in  an  overdetermined  system. 
The  motivation  to  this  solution  was  accuracy  improvement  and  the  possibility  of 
space  applications  since  an  initial  position-clock  bias  guess  was  not  required.  Ban¬ 
croft’s  solution  involves  solving  a  quadratic  equation,  where  each  of  the  two  roots 
leads  to  a  potential  solution,  one  of  which  does  not  satisfy  the  pseudorange  equa¬ 
tions  and  can  be  readily  eliminated.  This  solution  had  a  great  deal  of  merit  and 
motivated  a  number  of  papers  in  the  years  that  followed.  Driven  by  accuracy  and 
computational  issues  including  lower  dimensionality  and  speed,  Lloyd  Krause  [14] 
formulated  a  direct  solution  to  the  GPS  pseudorange  equation  of  the  exactly  deter¬ 
mined  system  based  on  difference  linearization.  By  differencing  the  satellite  position 
vectors,  a  new  basis  is  formed  by  using  any  two  adjacent  difference  vectors,  forming 
a  measurement  plane,  and  a  vector  orthogonal  to  the  plane.  The  four  nonlinear 
pseudorange  equations  expressed  in  the  new  basis  are  reduced  to  three  linear  equa¬ 
tions  that  are  independent  of  the  user  clock  bias  and  are  used  to  solve  for  the  user 
position  directly.  A  quadratic  auxiliary  equation  is  then  formed  to  solve  for  the  user 
clock  bias.  Krause’s  paper  demonstrated  a  brilliant  approach  by  which  differencing 
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is  used  to  linearize  quadratic  equations  and  remove  dependence  on  variables.  A  sim¬ 
ilar  approach  will  be  used  in  the  development  of  the  closed-form  algorithm  for  this 
thesis  research. 

Abel  and  Chalfee  [1]  demonstrated  that  in  both  closed-form  solutions  presented 
by  Bancroft  [2]and  by  Krause  [14],  a  position  fix  may  not  exist  and  if  it  does  exist,  it 
may  not  be  unique.  Abel  and  Chaffee’s  paper  concluded  that,  in  order  to  guarantee 
a  unique  position  fix,  an  overdetermined  system  using  at  least  five  satellites  must  be 
considered.  In  a  subsequent  paper  [6],  Abel  and  Chaffee  suggest  that  in  a  pseudo¬ 
range  system  such  as  GPS,  the  geometry  is  hyperbolic,  unlike  the  spherical  geometry 
of  a  ranging  system.  In  a  ranging  system,  ranges  are  measured  directly,  unlike  the 
case  for  GPS  in  which  the  pseudoranges  include  the  unknown  receiver  clock  bias. 
The  solution  to  range  equations  is  obtained  geometrically  through  the  intersection 
of  spheres,  but  this  method  does  not  generalize  to  pseudorange  equations  because 
of  the  unknown  bias  in  each  pseudorange;  hence,  it  is  not  possible  to  determine  the 
spheres.  In  view  of  the  fact  that  the  pseudoranges  are  not  only  corrupted  by  an 
unknown  clock  bias  but  also  by  some  measurement  noise,  caution  must  be  taken  in 
dealing  with  the  pseudorange  equations  when  it  comes  to  the  use  of  solutions  based 
on  spherical  geometry. 

2.2.4  Stochastic  Modeling.  Pseudorange  noise  that  corrupts  the  pseudo¬ 
range  measurements  is  caused  by  the  residual  errors  discussed  earlier.  In  order  to 
model  the  GPS  pseudorange  equations  statistically,  tremendous  effort  would  have  to 
be  dedicated  towards  the  development  of  reliable  noise  models.  This  noise  is  actually 
the  manifestation  of  receiver  noise  and  residuals  of  various  measurement  errors  that 
remain  unmodelled  and  uncompensated.  The  major  contributors  to  pseudorange 
noise  that  warrant  consideration  will  be  discussed. 
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Although  Gaussian-like,  receiver  noise  is  better  modeled  by  a  longer  tailed 
mixture  of  Gaussian  distributions  that  can  be  expressed  as: 

F(x)  =  (1  -  e)$(a;)  + 

where  $  is  a  Gaussian  distribution  and  the  parameter  e  is  generally  between  0.01  and 
0.1  [7].  In  a  multichannel  receiver,  receiver  noise  can  be  considered  to  be  uncorrelated 
across  satellites.  Consideration  must  also  be  given  to  the  receiver  clock  bias  which 
tends  to  correlate  the  pseudorange  measurements.  The  uncompensated  residuals  of 
tropospheric  and  ionospheric  errors  may  have  nonzero  means  that  also  add  to  the 
modeling  difficulties.  Noise  on  the  satellite  position  from  the  ephemeris  data  will 
also  have  some  effect  on  pseudorange  noise,  but  it  is  likely  to  be  non-Gaussian  and 
have  a  nonzero  mean. 

Given  the  number  of  contributing  factors  to  pseudorange  noise  and  our  lack 
of  knowledge  of  their  characteristics,  it  is  reasonable  to  propose  that  the  overall 
pseudorange  noise  will  have  a  zero-mean  Gaussian  distribution  by  invoking  the  Cen¬ 
tral  Limit  Theorem.  The  Central  Limit  Theorem  states  that  the  sum  of  many  inde¬ 
pendent  random  variables,  regardless  of  their  distribution,  will  approach  a  Gaussian 
distribution  [18].  The  Gaussian  pseudorange  noise  will  not  be  white  due  to  the 
correlated  nature  of  the  encompassed  errors  and  noise.  This  concern  is  alleviated 
since  there  is  no  requirement  for  the  pseudorange  measurements  to  be  uncorrelated 
in  time  because  the  positioning  problem  will  be  treated  as  a  static  estimation  prob¬ 
lem,  where  each  snapshot  in  time  is  treated  as  a  new  static  estimation  problem. 
On  the  other  hand,  it  is  desirable  to  the  development  of  the  stochastic  estimation 
that  pseudorange  noise  be  uncorrelated  across  satellites.  A  solid  argument  for  this  is 
not  available,  but  this  will  not  hinder  the  development  of  the  stochastic  closed-form 
solution  to  the  GPS  pseudorange  equations  in  this  thesis  since  the  pseudorange  mea¬ 
surements  will  be  differenced,  thereby  eliminating  some  of  the  effects  of  correlated 
noise.  The  uncorrelated  noise  after  differencing  can  only  be  justified  if  the  effects 
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of  Selective  Availability  (SA)  that  are  highly  time  correlated  in  nature  are  not  con¬ 
sidered  [30].  The  choice  to  overlook  the  effects  of  SA  in  this  thesis  is  not  overly 
restrictive  since  authorized  military  users  of  GPS  are  not  subject  to  the  effects  of 
SA  and  it  is  believed  that  the  use  of  SA  will  be  abolished  in  the  near  future. 

The  assumption  that  the  noise  across  pseudoranges  can  be  modeled  as  an 
independent  zero-mean,  variance  random  variable  is  used  by  Dailey  and  Bell  [8] 
without  any  solid  justification.  The  assumption  is  used  to  derive  statistics  for  the 
pseudorange  equations  position  solution  errors  similar  to  what  is  being  proposed  for 
this  thesis  research;  however,  their  approach  does  not  consider  a  closed-form  solution. 

2.2.5  The  Stochastic  Estimation  Problem.  The  components  of  an  estima¬ 
tion  problem  are  the  variables  to  be  estimated,  the  measurements,  and  a  mathemat¬ 
ical  model  describing  relationship  between  the  measurements  and  the  variables  to 
be  estimated  [18].  Given  the  lack  of  dynamics  in  the  GPS  pseudorange  equations 
at  any  given  time  instant,  a  static  estimation  problem  can  be  formulated  from  the 
stochastically  modeled  pseudorange  equations.  The  variables  to  be  estimated  and 
the  measurement  noise  on  the  pseudoranges  can  be  represented  by  random  variables. 
The  stochastic  estimation  process  will  not  only  provide  an  optimal  estimate  of  the 
unknown  variables,  the  user  coordinates  and  user  clock  bias;  but,  will  also  provide 
the  estimation  error  covariance.  This  is  the  most  significant  motivation  to  pursuing 
a  stochastic  approach  to  solving  the  GPS  pseudorange  equations.  However,  the  error 
covariance  accuracy  will  be  limited  by  the  quality  of  the  measurement  noise  statis¬ 
tics.  In  a  stochastic  estimation  problem,  this  emphasizes  the  requirement  for  good 
noise  models.  In  this  research  the  measurement  noise  is  modeled  as  a  zero-mean, 
(T^  variance  Gaussian  noise,  which  is  believed  to  be  adequate.  Due  to  the  lack  of 
knowledge  of  the  variance  a'^  and  the  fact  that  it  is  dependent  on  receiver  design, 
location,  orientation  and  time  of  day,  attempting  to  model  the  noise  variance  is  not 
a  viable  option.  The  approach  that  is  proposed  in  this  research  is  to  determine  the 
variance,  (cr^),  of  the  Gaussian  measurement  noise,  by  using  the  return  difference  or 
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measurement  residual.  It  is  suggested  that,  if  sufficient  measurements  are  available, 
reliable  data-driven  noise  statistics  can  be  extracted  by  evoking  ergodic  approxitions. 
The  question  at  hand  is  determining  how  many  measurements  are  sufficient  and  if 
that  is  a  reasonable  number  in  the  context  of  the  NAVSTAR  GPS  constellation. 

2.2.6  Conclusion  of  Literature  Review.  The  literature  review  supports 
the  proposed  approach  for  this  thesis  research.  There  are  problems  associated  with 
the  currently  used  iterative  approach  to  GPS  positioning  and  there  is  potential  for 
improvement  with  an  exact  closed-form  solution.  With  the  exception  of  one  recent 
paper  [8],  the  pseudorange  equations  have  generally  been  treated  as  a  determinis¬ 
tic  set  of  equations.  The  lack  of  effort  in  the  area  of  stochastic  modeling  applied 
to  GPS  pseudorange  equations  is  evident  from  the  lack  of  literature  on  the  sub¬ 
ject.  Previous  works  on  closed-form  solutions  for  the  GPS  pseudorange  equations 
did  not  make  use  of  the  pseudorange  measurements  from  all  in-view  satellites;  the 
derivations  considered  the  exactly  determined  case  using  only  four  of  the  available 
pseudorange  measurements  to  obtain  positioning  solutions.  Based  on  the  litera¬ 
ture  review,  the  development  and  evaluation  of  a  closed-form  solution  to  the  GPS 
pseudorange  equations  using  stochastic  modeling  and  estimation,  and  making  use  of 
all  in-view  satellites,  is  warranted. 
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III.  Mathematical  Derivations 

This  chapter  presents  a  thorough  development  of  the  closed-form  mathematical 
solution  to  the  GPS  pseudorange  equations.  The  development  of  supplementary  al¬ 
gorithms  that  make  lise  of  a  nonlinear  measurement  equation  and  that  were  deemed 
essential  to  enhance  the  performance  of  the  closed-form  algorithm  is  also  addressed. 
The  purpose  of  these  enhanced  algorithms  is  to  produce  an  improvement  to  the  solu¬ 
tion  achieved  by  the  closed-form  algorithm;  they  are  not  capable  of  providing,  nor  are 
they  intended  to  provide,  stand-alone  solutions.  The  two  stage  algorithms  that  are 
presented  are  the  Maximum  Likelihood  approach  and  the  Kalman  Update  approach. 
The  background  theory  and  the  operating  concept  behind  each  approach  is  discussed 
and  the  complete  mathematical  derivations  are  presented.  The  development  of  an 
alternate  closed-form  algorithm  is  also  presented. 

3.1  Stochastic  Closed-Form  Solution 

The  mathematical  derivation  of  the  stochastic  closed-form  solution  is  developed 
in  four  parts.  The  first  part  presents  the  algebraic  manipulations  to  transform  the 
stochastic  GPS  pseudorange  equation  as  shown  in  Equation  (1.1)  into  the  desired 
matrix  linear  regression  form  as  shown  in  Equation  (1.2).  The  second  part  involves 
the  derivation  of  statistics  for  the  equation  error  in  the  linear  regression.  The  third 
part  presents  the  development  of  the  static  stochastic  estimator  based  on  a  minimum 
variance  estimate  that  will  provide  an  estimate  of  the  user  position  coordinates 
and  user  clock  bias.  The  final  part  presents  the  derivation  of  the  estimation  error 
covariance  matrix. 

3.1.1  Basic  Concepts.  Prior  to  initiating  the  mathematical  derivation,  it 
is  necessary  to  present  some  basic  concepts  and  notation  that  will  be  used  in  the 
sequel; 
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•  X  N{fj,,cr'^)  is  the  notation  used  to  represent  a  random  variable  (a;)  that 
has  a  Gaussian  (Normal)  probability  distribution  function  with  mean  {n)  and 
variance 

•  E  is  used  to  represent  the  expectation  operator.  The  expectation  of  a  random 
variable  y  is  given  by  E{y}  =  pfy{p)dp,  where  fy{p)  is  the  density  of  y  [18]. 
This  also  defines  the  mean  (^),  the  first  moment  of  the  random  variable. 

•  Expectation  is  a  linear  operation;  therefore,  for  any  two  random  variables  x 
and  y,  then  E{x  +  j/}  =  E{x}  +  E{y}. 

•  If  two  random  variables  x  and  y  are  uncorrelated,  then  E{xy}  =  E{x}E{y}. 

•  The  variance  of  the  element  of  a  random  vector  x  can  be  expressed  as. 
Pa  =  E{{xi  —  PiY}i  which  make  up  the  diagonal  elements  of  the  covariance 
matrix.  The  off-diagonal  elements  Pij  are  zero  if  the  and  elements  of 
the  random  vector  x  are  uncorrelated. 

•  For  a  random  variable  x  ~  N{0,  <t^),  the  moments  of  x  are  expressed  as  E{x’^} 

for  k  =  (1, 2,  •  •  • ,  oo).  An  odd  k  denotes  an  odd  moment  and  an  even  k  denotes 
an  even  moment.  All  odd  moments  are  zero  and  even  moments  are  given  by 
{k  -  [23]. 

8.1.2  Linear  Regression.  The  derivation  of  the  closed-form  solution  to  the 
pseudorange  equations  begins  with  Equation  (1.1).  The  general  form  solution  with 
n  >  5  satellites  is  presented. 

Ri  -  \J (Ua,  -  XiY  -1-  {Uy  -  yiY  {u^  -  ZiY  -\-b  +  Wi, 

represents  the  pseudorange  equation  of  n  equations  where  Wi  is  Gaussian  mea¬ 
surement  noise  with  statistics  Wi  ~  N{G,a‘^)  and  E{wiWj}  =  0  for  (z  ^  jf).  The 
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pseudorange  equation  can  also  be  rewritten  as: 

{u^  -  Xif  +  {Uy  -  yif  +  {Uz  -  Zif  =  {Ri-b-  Wif 

By  expanding  both  sides  of  the  above  equation  and  rearranging  the  terms  such  that 
only  the  known  terms  and  the  noise  terms  are  on  the  left  hand  side,  the  equation 
becomes: 

"I"  Uy  d"  U^  b  ‘IXiUx  ^Vi'^y  ^ZjUn  QiRib  — 

i?-  -  X-  -  Ui  -  zf  -  2RiWi  +  2bwi  +  wf  (3.1) 

It  is  noted  that  the  first  four  terms  in  Equation  (3.1)  are  simply  the  unknown  vari¬ 
ables  squared  and  that  they  are  common  to  all  n  equations.  This  presents  an  oppor¬ 
tunity  of  eliminating  the  nonlinear  terms  by  differencing;  hence,  the  equation  is 
subtracted  from  the  remaining  n  —  1  equations.  The  resulting  n  —  1  equations  are 
linear  in  the  unknown  variables  and  can  be  expressed  as: 


(Xn  Xj)'Uj;  -(-  (t/n  T  (•^n  Z^U^  T  {^Ri  Rn)b 


i(i2-  -Rl  +  xl-xj  +  yl- yf  ■\-zl-z1)  +  R^Wn - RiWi  +  bwi - bwn  +  ^(u^f - wl)  (3.2) 

As  a  by-product  of  the  preceding  operation,  the  nonlinear  pseudorange  equation 


remains. 

Rn  =  \J{Ux-  XnY  +  {Uy  -  t/n)^  +  (Wz  -  Z^Y  +  b Wn  (3.3) 


The  equation  will  remain  unused  for  this  section  of  the  derivation  but  will  be  used 
subsequently  as  an  auxiliary  equation  for  developing  a  maximum  likelihood  solution 
and  for  use  in  the  Kalman  update  solution. 

To  simplify  Equation  (3.2)  further,  an  average  pseudorange  expression  is  in¬ 
troduced  by  defining: 

R=-j2Ri  (3-4) 
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The  following  approximation  is  also  introduced  as  a  first  attempt  at  reducing 
Equation  (3.2)  by  simplifying  the  noise  terms  RnWn  —  RiWi’. 

Ri^Rn^R  (3.5) 

In  the  GPS  constellation  it  is  not  possible  to  produce  a  scenario  with  more  than 
five  satellites  in  which  all  pseudoranges  are  the  same;  therefore,  the  approximation 
in  Equation  (3.5)  introduces  a  source  of  error  whereby  some  geometrical  scenarios 
will  be  more  prone  to  errors  than  others.  The  impact  of  introducing  the  approxi¬ 
mation  is  not  known  but  it  will  undoubtedly  affect  the  accuracy  of  the  closed-form 
solution  being  developed.  If  the  effect  is  found  to  be  significant  through  experimen¬ 
tal  analysis,  this  approximation  must  be  revisited.  To  minimize  the  impact  of  the 
approximation  in  Equation  (3.5),  the  measured  pseudoranges  will  be  ordered  such 
that  the  pseudorange  closest  to  R  will  be  used  as  the  measurement.  By  intro¬ 
ducing  the  approximation  in  Equation  (3.5)  into  the  noise  terms  of  Equation  (3.2), 
the  equation  becomes: 


{Xffi  xfjVjx  -j-  (j/n  d~  (^n  d”  ^Ri  Rjijl^  ^ 


^(^i  +  +  yf  +  4  -  ^i)  -  -  ^n)  +  -  Wn)  +  -  wl) 


From  this  point  on  in  the  derivation,  the  exactness  of  the  solution  has  been  lost  due 
to  the  approximation  that  has  been  introduced.  The  equation  is  further  reduced  to 
obtain: 


(a:„  -  Xi)Ux  -I-  (j/n  -  yi)Uy  -I-  {Zn  -  Zi)u^  -t-  {Ri  -  Rn)h  W 
i(i??  -  Rl-^xl-x]-\-yl-  yf  +  zl-  zf)  +  (6  -  R){wi  -  w^)  +  ^{wf  -  wf) 


The  following  noise  terms  are  redefined  to  allow  for  further  simplification  of  the 
equation: 
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•  {Wi  -  Wn)  =4»  Vi 

The  equation  can  then  be  rewritten  as: 

(^n  T  (2/71  Vi^'^y  ”b  (-^n  T  (-^z 

+  2/n  -  4  -  ^i)  +  (^  - 

By  redefining  the  following  terms: 

•  (b  —  R)vi  +  Wi=^  K-; 

•  +  +yl- Vi  +4-  A)  ^ 

the  equation  reduces  to  a  linear  regression  form  of  (n  — 1)  equations  in  four  unknowns 
and  is  expressed  as: 

(^n  ^i')Vx  "}“  2/z)^2/  ”1“  (^Ri  Rn)b  ^  T  kj  (3.6) 

The  linear  regression  in  Equation  (3.6)  can  be  compactly  written  in  matrix  notation 
form  as: 

Z  =  HX  +  V,  (3.7) 

where 

•  .Z  is  the  (n  —  l)-dimensional  measurement  vector  given  by: 

Z2 

z=  : 

Zn-l 
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•  if  is  the  (w  —  1)  X  4  regressor  matrix  given  by: 

^1  Vn  Vl  ^1  i^l 

^2  Vn  2/2  ^2  -^2  i^n 

H=  \  \  \  \ 

3Jji  Xji—l  yn  Vn—l  ^n— 1  ^n—1  i^n 

•  X  is  the  vector  of  unknowns  given  by: 

y>x 

X  = 

Uz 

h 

•  y  is  the  {n  —  1)  equation  error  vector  given  by: 

y 
y 

V= 

_  K-i 

S.l.S  Noise  Statistics.  Using  the  linear  regression  obtained  in  Section  3.1.2 
as  shown  in  Equation  (3.7),  it  is  possible  to  obtain  an  estimate  of  X  if  the  statistics 
of  the  noise  vector  V  are  given.  The  noise  vector  statistics  are  not  known;  however, 
the  composition  of  the  noise  vector  elements  is  known.  The  statistics  of  V  must  be 
derived  from  the  known  statistics  of  the  pseudorange  measurements  noise,  w,-. 

For  the  statistical  derivations  that  follow,  when  the  subscript  i  and  the  sub¬ 
script  j  are  used  together  it  is  implied  that  (i  ^  j). 
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Statistics  for  w  where  (i,j  =  l,2,---,n): 


E{wi}  =  0 

E{w^}  = 

E{wiWj}  = 

Statistics  for  v  and  w  where  {i,j  =  1, 2,  •  •  • ,  n  —  1)  : 

Vi  =  {Wi  -  Wn) 

Wi  = 

E{Vi}  =  E{{Wi-Wn)} 

=  E{W{}  -  E{Wn} 

=  0-0  =  0 


E{vf}  =  E{{Wi-WnY} 

=  E{wi}  4-  E{wl}  -  2E{wiWn} 

=  (T  V  (7^  -  0  = 


E{viVj}  =  E{{wi  -  Wn)iwj  -  u;„)} 

=  E{wiWj}  -  E{wiWn}  —  E{wjWn}  +  E{wl} 
=  0  —  0  —  0  +  <7^  =  (T^  (i  j) 
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E{w,}  = 

=  \(E{>^n  -  E{«i]) 

=  \(a^-a^)  =  0 


Eiwl)  =  E{{\(w^-wl))^} 

=  ^{E{wt]  +  E{wt,]-2E{w^u,l}) 

=  7(3<7^  +  3(7'*  —  2(7*)  =  a* 

4 

E{wiWj}  = 

=  -  <7^^  -  (7^)  =  i(7^  (i  ^  i) 

=  E{{wi  -  Wn)^iw^  -  wl)} 

=  \{E{w'^i}  +  E{wl}  -  E{wiwl}  -  Eiw'lwn}) 
=  ^(0  +  0-0-0)  =  0 

E{viWj}  =  E{{wi-wn)]^{w]-wl)] 

=  ^{E{wiw]]  +  E{wl]  -  Eiwiwl)  -  E{w]wn}) 

=  i(0  +  0-0-0)  =  0(i^i) 


3-8 


Statistics  for  V  where  (i,  j  =  1, 2,  •  •  • ,  n  —  1)  : 


Vi  =  {b-  R)vi  +  Wi 


E{Vi}  =  E{{b-R)vi  +  Wi} 

=  (b  —  R')E{vi}  +  E{wi} 

=  0+0=0 


E{V^}  =  E{{{b~R)vi  +  Wi){{b-R)vi  +  Wi)} 

=  (6  -  RYE{vf}  +  2(6  -  R)E{viWi}  +  E{w]} 

=  (6  -  RY(2a^)  +  2(6  -  ^)(0)  +  <7^ 

=  2(6  -  Rf(x^  +  (7^  (3.8) 


E{ViVj}  =  E{i{b-R)vi  +  Wi){{b-R)vj  +  Wj)} 

=  (6  —  .R)^  +  (6  —  R)E{viw;j}  +  (6  —  R)E{vj‘Wi}  +  E{wiWj} 

=  (6  -  .R)^(<7^)  +  (6  -  ^)(0)  +  (6  -  ^)(0)  + 

=  (6  -  (i  ^  i)  (3.9) 

From  the  results  obtained  in  Equation  (3.8)  and  Equation  (3.9)  it  can  be  seen 
that  the  covariance  of  the  noise  vector  V  has  diagonal  elements  that  are  twice  the 
off-diagonal  elements.  Indeed,  the  covariance  matrix  of  the  equation  error  vector,  V, 
can  be  expressed  as: 

Pv  =  (r\^(T^  +  {b-RY)Rn,  (3.10) 
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where  is  a  (n  —  1)  x  (n  —  1)  matrix  defined  as  follows: 


2 

1 

1 

...  1 

1 

2 

1 

...  1 

Rn  = 

1 

1 

2 

1 

1 

2 

(3.11) 


An  interesting  characteristic  of  the  matrix  is  that  its  determinant  is  always  equal 
to  n  which  allows  it  to  be  easily  inverted  analytically.  The  inverse  of  i2„  is  required 
to  find  the  inverse  covariance  matrix  that  is  given  by: 


p-i  _ _  _ R-i 

^  (T2(i(T2  +  (6-5)2)  ’ 


and  where  is  the  matrix 

inverse 

of  Rn 

explicitly 

given  by: 

n  —  1 

-1 

-1  • 

..  _i 

-1 

n  —  1 

-1  • 

..  _i 

-1 

-1 

n  —  1 

: 

n 

; 

: 

; 

-1 

-1 

•  «  . 

•  •  n  —  1 

(3.12) 


(3.13) 


3.14  Minimum  Variance  Estimate  Solution.  Using  the  linear  regression 
from  Equation  (3.7)  as  a  starting  point,  the  aim  is  to  obtain  an  estimate  X  of  the 
vector  of  unknown  parameters  X.  The  X  that  minimizes  the  estimation  error  as 
weighted  by  the  inverse  covariance  of  the  noise  must  be  obtained.  Recognizing  that 
the  estimation  error,  also  known  as  the  return  difference  or  measurement  residual, 
is  (Z  —  HX)j  the  estimation  problem  can  be  formulated  as  follows: 

mm^[(Z  -  HXfPv\Z  -  HX)]  (3.14) 
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Equation  (3.14)  can  be  expanded  to  obtain: 


Z^Pv^Z  +  X^H^Py^HX  -  2X^H^Py^Z  (3.15) 

Since  a  minimization  over  X  is  needed,  Equation  (3.15)  is  differentiated  with  respect 
to  X  and  set  equal  to  zero  yielding  the  following  expression: 

0  +  2H^Pv^HX  -  2H^Py^Z  =  0 

Rearranging  the  expression  and  solving  for  X  produces  the  desired  solution: 

Jl  =  {H^Py^Hy'^H^Py^Z  (3.16) 

In  order  to  demonstrate  that  the  stationary  point  solution  in  Equation  (3.16)  is 
indeed  a  minimum,  the  hessian  matrix  of  Equation  (3.15)  must  be  verified.  The  re¬ 
sulting  hessian  matrix  is  (H^ Py^H),  which  by  definition  is  always  positive  definite, 
providing  the  necessary  and  sufficient  conditions  for  minimization.  Furthermore, 
since  the  existance  of  the  hessian  inverse  in  Equation  (3.16)  is  guaranteed,  the  exis- 
tance  of  a  solution  is  also  guaranteed. 

Equation  (3.16)  is  a  closed-form  solution  to  the  GPS  pseudorange  equations. 
The  solution  is  not  dependent  on  cr,  the  pseudorange  measurement  noise  standard  de¬ 
viation.  To  simplify  the  solution  for  implementation,  it  is  noted  that  Equation  (3.10) 
shows  the  noise  covariance  Py  as  simply  premultiplied  by  a  scalar  quantity.  In 
Equation  (3.16)  the  scalar  premultiplier  of  Py  will  cancel  out;  therefore,  the  esti¬ 
mation  solution  shown  in  Equation  (3.16)  can  be  rewritten  in  an  equivalent  form 
as: 

1  =  {H^RZ^Hy^H^PyZ  (3.17) 

Equation  (3.17)  is  used  for  coding  the  experimental  Matlab  algorithm.  It 
must  be  noted  that  there  are  no  big  matrix  inversions  associated  with  this  solution 
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since  R~^  has  been  determined  analytically  and  can  be  coded  directly  into  the  algo¬ 
rithm.  The  only  inversion  that  needs  to  be  performed  is  that  of  the  (4  x  4)  matrix, 
{H^ R~^H),  which  can  be  hardwired  into  the  receiver’s  algorithm. 

S.1.5  Estimate  Error  Covariance.  It  follows  from  Equations  (3.14)  and 
(3.16)  that  the  covariance  of  the  estimate  error  is  given  by: 

:px  =  E{{x  -  j()(x  -  kf] 

=  +  (i  -  (3.18) 

Unlike  the  solution  estimate,  the  covariance  Px  is  dependent  on  (t;  hence,  a  must  be 
known  or  estimated,  in  order  to  compute  the  error  covariance.  If  a  sufficiently  large 
number  of  GPS  pseudorange  measurements  are  available,  the  following  approach  of 
using  an  ergodic  assumption,  whereby  an  ensemble  average  is  approximated  by  a 
time  average,  can  attempted: 

e{Av}  = 

=  (n  -  1)B{K"} 

=  (n  —  1)(2(6  —  -b  u"*)  (3.19) 

Using  the  return  difference  from  the  measurement  data  and  using  an  ergodic  ap¬ 
proximation,  an  expression  equivalent  to  Equation  (3.19)  can  also  be  obtained  as 
follows: 

E{V^V]k{Z -HXf{Z -Hk)  (3.20) 

By  equating  Equations  (3.19)  and  (3.20)  and  re-arranging  the  terms,  a  quadratic 
equation  in  <7^  is  obtained: 

a'^^2{h-Rfa'^-^—{Z-Hkf{Z-HX)  =  ti 

Tt  x 
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Solving  the  quadratic  equation  yields  the  following  data  driven  estimate  of  cr^: 


=  -{R  -  hf  +  w  -  6)"  +  -  Hkyiz  -  nk)  (3.21) 

y  Tt  1. 

If  Equations  (3.19)  and  (3.20)  are  accepted  as  being  equivalent  through  evoking 
ergodicity,  the  expression  for  the  error  covariance  matrix  given  in  Equation  (3.18) 
can  be  rewritten  in  terms  of  the  measurement  data  driven  return  difference  as  follows: 

Px  =  -  Hkf(Z  -  (3-22) 

It  must  be  recognized  that  the  validity  of  Equation  (3.22)  is  dependent  on  the 
validity  of  the  ergodic  approximation  introduced  at  Equation  (3.20);  hence,  a  yet-to- 
be-determined  minimum  satellite  availability  (n^m)  must  be  present  to  validate  this 
solution  for  the  estimation  error  covariance.  An  acceptable  value  for  (rimin)  must  be 
determined  experimentally. 

3.2  Maximum  Likelihood  Estimate 

This  section  presents  the  background  theory  behind  the  maximum  likelihood 
solution  and  a  potential  estimation  algorithm  is  developed.  It  is  important  to  note 
that  the  solution  that  is  presented  in  this  section  is  not  a  stand-alone  solution  to 
the  GPS  pseudorange  equations;  it  is  a  supplementary  process  that  improves  on  the 
closed-form  solution  presented  in  Section  3.1. 

3.2.1  Maximum  Likelihood  Concept.  The  concept  behind  the  maximum 
likelihood  solution  approach  is  simple  and  ideally  suited  as  a  natural  continuation  to 
the  closed-form  solution.  Given  a  known  probability  density  function  of  a  random 
vector  and  some  equality  constraint(s)  on  the  random  vector,  it  is  possible  using  op¬ 
timization  techniques  to  find  the  maximum  value  of  the  density  function  subject  to 
the  constraint  (s).  The  optimized  maximum  value  of  the  probability  density  function 
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corresponds  to  the  most  likely  realizations  of  the  random  vector.  Hence  the  solution 
approach  will  be  referred  to  as  Maximum  Likelihood.  The  concept  is  clearly  demon¬ 
strated  in  the  two  dimensional  illustration  in  Figure  3.1.  Figure  3.1  shows  a  Gaussian 
probability  distribution  function  in  two  variables  subject  to  a  linear  constraint  on 
the  two  variables.  Given  the  constraint,  the  values  for  the  two  random  variables 
that  maximize  the  probability  function  is  readily  observed  as  the  intersection  of  the 
constraint  with  the  highest  constant  probability  contour  line. 


-4  -4 

Figure  3.1  Maximum  Likelihood  of  Gaussian  Probability  Distribution  Subject  to  a 
Linear  Constraint 


The  closed-form  solution  to  the  GPS  pseudorange  equations  developed  in  Sec- 
tion  3.1  of  this  thesis  proposed  that  the  four  dimensional  random  vector  X  has  a 
Gaussian  probability  density  function.  The  closed-form  solution  provided  the  mean 
X  and  covariance  Px  which  can  be  represented  as: 

X  -  iv(i,  Px), 
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or  in  expanded  form  as: 


Ux,  Uy,  Uz,  bf  ~  iV([Uj;,  Uj,,  hf,  Px). 


A  Gaussian  probability  density  function  is  fully  described  by  its  first  two  mo¬ 
ments;  hence,  although  it  is  impossible  to  visualize  the  four  dimensional  density 
function  in  our  three  dimensional  world,  it  can  be  represented  mathematically  along 
with  an  equality  constraint  to  formulate  an  optimization  problem. 


3.2.2  Mathematical  Development.  The  Gaussian  probability  density  func¬ 
tion  of  a  A:-dimensional  random  vector  (T)  whose  realizations  are  represented  by  the 
vector  7  is  described  by  [18]: 


fvil)  = 


1 

(27r)t|Pl? 


(3.23) 


where  P  is  the  {k  x  k)  covariance  matrix  and  m  is  the  {k  x  1)  mean  vector.  It  follows 
from  Equation  (3.23)  that,  in  the  context  of  the  GPS  solution  estimate,  the  Gaussian 
probability  density  function  for  the  random  vector  (X)  can  be  expressed  as: 


1 

(27r)2|P;,|^ 


(3.24) 


The  goal  in  this  section  is  to  find  the  GPS  solution  (X)  that  maximizes  the 
probability  density  function  f{X)  in  Equation  (3.24).  In  order  to  maximize  the 
function  f{X),  {X)  must  be  selected  such  that  the  following  part  of  the  exponential’s 
exponent  is  minimized: 

[x-xfp^'ix-k] 
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The  exponent  to  be  minimized  can  be  rewritten  in  expanded  form  as: 

r  T  ^  r  n 


Ujj  Ux 

Uy  Uy 

p-l 

Uy  Uy 

(3.25) 

U^ 

Uz  ^2: 

b-b 

A 

b-b 

At  this  point  it  is  necessary  to  introduce  the  pseudorange  measurement 
equation  described  in  Equation  (3.3).  The  exponent  term  in  Equation  (3.25)  must 
be  minimized  subject  to  the  constraint  imposed  by  Equation  (3.3).  Equation  (3.3) 
cannot  be  used  as  a  constraint  in  its  present  form  due  to  the  noise  term  it  contains. 
The  zero-mean  noise  term,  in  the  pseudorange  measurement  equation  is 
neglected  to  obtain  a  simple  equality  constraint.  A  minimization  problem  can  now 
be  formulated  ais  follows: 


subject  to: 


n^i'>^(Ux,Uy,Uz,b) 


1 

■ 

Ux  Ux 

Ux  Ux 

Uy  Uy 

p-l 

1 

Uz  Uz 

Uz  Uz 

1 

1 

_ 1 

b-b 

Rn  —  \l  (^a 


“I”  (Uy 


VnY  +  {Uz  -  ZnY  +  b 


(3.26) 


(3.27) 


The  constrained  optimization  problem  specified  in  equations  3.26  and  3.27  can 
be  solved;  however,  it  is  desirable  to  use  the  equality  constraint  to  solve  explicitly  for 
one  of  the  optimization  parameters,  thereby  reducing  the  order  of  the  optimization 
problem  as  suggested  by  Maybeck  [19].  Equation  (3.27)  can  be  rearranged  to  solve 
for  the  user  clock  bias  (&): 

h=  Rn  -  ^ {Ux  -  Xnf  +  {Uy  -  VnY  +  {Uz  -  ZnY  (3.28) 
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Equation  (3.28)  should  be  showing  a  non-unique  solution  for  b  caused  by  the 
positive  and  negative  roots;  however,  the  positive  root  has  been  dropped  since  it 
leads  to  a  physically  impossible  solution.  Substituting  Equation  (3.28)  into  the  min¬ 
imization  problem  defined  in  Equation  (3.26)  will  reduce  the  optimization  problem 
to  an  unconstrained  minimization  in  three  variables  which  can  be  solved  with  greater 
ease. 


^^^{Ux,Uy,Uz) 


Uz  - 


T 


^x  '^x 

Uy  Uy 

Uz  —  Uz 

_  (Rn  -  ‘sjiv-x  -  XnY  +  {Uy  -  VnY  +  (Wz  “  ^n)^)  “  & 

After  Equation  (3.29)  is  solved  for  the  optimization  parameters  {uxiUy,Uz)  which 
correspond  to  the  Meiximum  Likelihood  GPS  user  position  estimate,  (ua,,  Uy,  Uz)  must 
be  substituted  into  Equation  (3.28)  to  obtain  the  Maximum  Likelihood  GPS  user 
clock  bias  estimate  (6). 

The  unconstrained  minimization  problem  in  Equation  (3.29)  can  be  solved  us¬ 
ing  a  number  of  different  optimization  approaches.  Solving  the  actual  optimization 
problem  is  beyond  the  scope  of  this  thesis;  however,  the  function  used  for  the  exper¬ 
imental  implementation  of  the  Maximum  Likelihood  solution  needs  to  be  discussed. 
The  experimental  algorithm  implemented  in  Matlab  [17]  produced  the  most  desirable 
results  using  the  function  fmins.  The  search  initialization  point  used  experimentally 
was  {uxi  Uy,  Uz)  corresponding  to  the  previous  best  GPS  position  solution  estimate 
provided  by  the  closed-form  solution  from  Section  3.1.  The  fmins  Matlab  function 


(3.29) 
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uses  the  simplex  search  method  of  Nelder  and  Mead  [17].  The  simplex  search  method 
does  not  use  any  derivative  or  gradient  information,  rather  it  involves  the  use  of  a 
simplex.  In  m-dimensional  space,  where  m  corresponds  to  the  order  of  the  optimiza¬ 
tion  problem,  the  simplex  is  characterized  by  (m  -)- 1)  distinct  vertices.  At  each  step 
in  the  search,  one  of  the  vertices  is  replaced  by  a  new  point  reducing  the  size  of  the 
simplex.  The  search  continues  until  the  size  of  the  simplex  is  less  than  the  stopping 
tolerance  [17].  The  fmins  function  consistently  produces  a  solution  to  the  minimiza¬ 
tion  problem  regardless  of  how  good  the  previous  estimate  used  for  initialization  was. 
The  fminu  function  which  uses  a  quasi-Newton  search  method  was  also  tried  in  the 
experimental  algorithm  but  it  was  not  as  robust  as  the  fmins  function.  The  fminu 
function  occasionally  failed  to  produce  a  minimum  when  the  search  initialization 
point  was  not  close  enough  to  the  true  user  position. 

The  Maximum  Likelihood  approach  described  makes  effective  use  of  the  stochas¬ 
tic  estimation  solution  produced  by  the  closed-form  algorithm  presented  in  Sec¬ 
tion  3.1  and  of  the  pseudorange  equation,  that  had  not  previously  been  used, 
to  provide  a  better  GPS  solution  estimate.  If  the  constraint  formed  from  the 
pseudorange  equation  were  truly  a  deterministic  equation  as  was  assumed  in  Equa¬ 
tion  (3.27),  the  Maximum  Likelihood  approach  would  produce  the  true  user  position 
and  user  clock  bieis  as  its  GPS  solution  estimate.  Unfortunately,  the  pseudo¬ 
range  equation  is  noise  corrupted,  as  are  the  other  {n  —  1)  pseudorange  equations; 
therefore,  the  Maximum  Likelihood  estimate  is  prone  to  the  effects  of  the  noise  on 
the  pseudorange  measurement.  Perhaps  the  largest  disadvantage  of  taking  a 
maximum  likelihood  approach  is  that  there  is  no  simple  way  of  determining  the  error 
covariance  of  the  the  solution  estimate  it  produces.  The  capability  to  estimate  the 
error  covariance  is  traded  off  to  obtain  an  improved  GPS  solution  estimate. 
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3.3  Kalman  Update  Solution 

This  section  presents  the  development  of  the  Kalman  update  GPS  position 
determination  algorithm.  As  was  the  case  for  the  maximum  likelihood  algorithm 
developed  in  Section  3.2,  the  Kalman  update  solution  is  not  a  stand-alone  solution 
to  the  GPS  pseudorange  equations;  it  is  a  supplementary  process  that  improves  on 
the  closed-form  solution  presented  in  Section  3.1.  The  concept  behind  the  Kalman 
Update  solution  is  discussed  followed  by  the  complete  mathematical  derivation  of  the 
solution  itself.  The  Kalman  update  algorithm  presented  in  this  section  is  an  alternate 
to  the  maximum  likelihood  algorithm  to  enhancing  the  GPS  solution  produced  by 
the  closed-form  algorithm  by  using  the  previously  unused  pseudorange  equation. 
However,  the  Kalman  update  algorithm  is  capable  of  obtaining  the  estimation  error 
covariance  which  is  not  readily  obtainable  using  the  maximum  likelihood  algorithm. 

3.3.1  Kalman  Update  Concept.  The  concept  behind  the  Kalman  update 
solution  approach  is  similar  to  that  of  a  conventional  Kalman  Filter.  The  closed- 
form  solution  in  Section  3.1  provides  a  GPS  solution  estimate  (X)  and  the  associated 
estimation  error  covariance  matrix  (Ux)-  Recalling  that  this  solution  was  produced 
without  making  use  of  the  pseudorange  equation  in  Equation  (3.3),  the 
pseudorange  equation  can  be  perceived  as  a  new  measurement  which  can  be  used  to 
update  the  previous  estimate  the  same  way  that  it  would  be  accomplished  during 
the  update  cycle  of  an  extended  Kalman  Filter.  The  approach  that  is  used  begins 
with  the  linearization  of  Equation  (3.3)  about  a  nominal  position  estimate.  The 
linearized  equation  is  then  manipulated  into  the  standard  linear  measurement  form 
as  described  in  [18],  and  used  to  update  the  estimate.  It  may  be  necessary  for  the 
process  to  continue  in  an  iterative  manner  until  convergence  within  a  predefined 
tolerance  is  achieved;  however,  it  is  expected  that  the  algorithm  will  converge  to  a 
solution  in  a  few  iterations. 
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The  Kalman  update  algorithm  that  is  presented  in  this  section  differs  from  the 
basic  Kalman  Filter  developed  by  Kalman  [12]  [13]  in  that  the  measurement  that  is 
used  to  update  the  previous  estimate  is  correlated  with  the  previous  estimate.  The 
conventional  Kalman  Filter  update  equation  does  not  allow  for  correlation  between 
the  new  measurement  and  the  previous  estimate;  hence,  a  Kalman-like  update  equa¬ 
tion  that  can  account  for  this  correlation  is  presented  for  use  in  the  experimental 
algorithm. 

S.S.2  Linearized  Measurement.  The  first  step  in  the  mathematical  develop¬ 
ment  of  the  Kalman  update  algorithm  is  to  linearize  Equation  (3.3)  about  a  nominal 
user  position  performing  a  Taylor  series  expansion  and  neglecting 

second  and  higher  order  terms.  The  linearized  equation  obtained  is  given  by 


By  defining  the  regressor  h  for  this  scalar  measurement  equation  as  follows: 


the  linearized  equation  can  be  rewritten  as 
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Rn  ^  h'^X  +  Wn 

(^a;o  ^n)^xo  4”  (^yo  UTi^'^yo  4”  (^^o  ^7i)^zo 

4"  (^t/o  Vn)^  4“  {}^zq 

4-  ^(l£a;o  -  +  (%o  ““  2/n)^  +  {^z^  “ 


where  X  is  the  vector  of  unknowns,  [ux-^Uy^Uz^hY •  The  goal  is  to  reduce  the  above 
equation  into  the  form  of  a  linear  measurement  model  described  by: 

Z,  =  h^it  +  u;,  (3.31) 


In  order  to  achieve  this  goal,  Zn  must  be  defined  as: 


Zn  —  Rn  {}^xq  4"  (^yo  J/n)^  4“  {}^zq 

^  (^370  ^n)'^XQ  4~  (^yo  i/n)^yo  4~ 

\l{}^xo  4*  (^yo  VnY  4“  {^'Zq  -^n) 

which  can  be  simplified  into  the  following: 

ry  r>  I  (^^0  ^n)^n  4"  (“^yo  2/n)2/n  4"  {"^zq  ^n)^n 

\  j - ^ — 

{y^XQ  4"  (’^yo  Vn)  4“  {^'^zq  -^n) 


(3.32) 


Now  that  the  pseudorange  measurement  equation  is  approximated  into  the 
appropriate  linear  measurement  model  form  defined  in  Equation  (3.31),  it  can  be 
used  to  update  the  solution  obtain  from  the  closed-form  algorithm  using  a  Kalman- 
like  update  approach.  Using  a  linear  measurment  model  simplifies  the  solution  by 
allowing  the  use  of  linear  Kalman  filtering  techniques,  as  apposed  to  using  an  Ex¬ 
tended  Kalman  filter  or  increasing  the  order  of  the  filter  to  accomodate  a  nonlinear 
measurement  equation.  Keeping  in  mind  that  Zn  is  actually  part  of  the  measure¬ 
ments  that  were  used  to  obtain  the  closed-form  solution  and  not  a  new  measurement 
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as  would  be  the  case  in  a  conventional  Kalman  Filter  application,  hence  the  new 
measurement  and  the  previous  estimate  are  correlated.  This  is  a  violation  to  the 
basic  assumptions  used  in  the  derivation  of  the  conventional  Kalman  Filter  update 
equations.  A  Kalman-like  update  equation  that  can  accommodate  correlation  be¬ 
tween  the  new  measurement  and  the  previous  estimate  needs  to  be  derived. 

S.S.S  Noise  Statistics.  In  order  to  derive  the  new  Kalman- like  update 
equation,  it  is  necessary  to  know  the  relationship  between  the  noise  in  the  new 
measurement  (i«„)  and  the  previous  estimate  being  the  solution  obtained  from  the 
closed-form  algorithm.  The  linear  regression  used  for  the  closed-form  algorithm  was 
defined  in  Equation  (3.7)  as, 

Z^HX  +  V, 

and  the  statistics  of  the  noise  vector  V  were  derived  in  Section  3.1.3.  The  closed- 

— f 

form  algorithm  produced  an  estimate  of  the  GPS  unknown  parameter,  (AT),  defined 
in  Equation  (3.17)  and  an  estimate  of  the  covariance  matrix  associated  with  estima¬ 
tion,  {Px)i  defined  in  Equation  (3.22).  Using  the  knowledge  of  the  estimated  GPS 
solution,  the  true  GPS  parameter  vector  can  be  defined  as, 

X=k  +  w,  (3.33) 

where  W  ~  A^(0,  P^)-  The  correlation  of  interest  between  and  W  can  be  defined 
as: 

p  =  E{Wwn}  =  E{wnW}  (3.34) 

— ^  »  • 

To  determine  the  relationship  between  W  and  V,  the  linear  regression  in  Equa¬ 
tion  (3.7)  is  multiplied  from  the  left  by  H^R~^  yielding  the  following  expression: 


The  expression  can  be  solved  for  X  to  obtain: 

X  =  -  {H^R-^H)-^H^RZ^V  (3.35) 

The  first  term  on  the  right  hand  side  of  Equation  (3.35)  is  recognized  from  Equa¬ 
tion  (3.17)  as  X]  therefore,  by  equating  Equation  (3.35)  and  Equation  (3.33)  an 
expression  for  W  in  terms  of  V  is  obtained: 


W  =  (H^RZ^H)-^H^R-^V 


(3.36) 


Next  the  relationship  between  V  and  Wn  is  determined  be  exploiting  the  noise 
statistics  derived  in  Section  3.1.3. 


E{ViWn}  =  E{{{b- R)Vi  +  Wi)Wn} 

=  {b-  R)E{viWn}  +  E{WiWn} 

1 


(6  -  R)E{{Wi  -  u;„)wn}  +  ^E{{Wi  -  wl)Wn} 

=  {b-  R)E{WiWn  -  wl}  +  ^E{WiWn  -  u^^)} 

=  (6-Ji)(0-<T")  +  i(0-0) 

=  {R-b)a^. 


(3.37) 


Equation  (3.37)  which  represents  the  variance  between  any  single  element  of  V  and 
Wn-  The  relationship  can  be  generalized  to  obtain  the  following  covariance  matrix: 


E{Vwn}  =  {R-  by 


(3.38) 


-I  {(n-l)xl) 
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Using  Equation  (3.36)  and  the  relationship  in  Equation  (3.38),  an  expression 
for  the  covariance  between  W  and  Wn  is  determined: 


=  {R- 


1 

1 


1 


((n-l)xl) 


(3.39) 


3.3.4  The  Augmented  Linear  Regression.  An  augmented  linear  regres¬ 
sion  can  be  formulated  by  combining  Equation  (3.33)  and  Equation  (3.31).  The 
augmented  linear  regression  is  expressed  as: 


Za  =  HaX  +  Va 


(3.40) 


where, 

•  Za  is  the  (5x1)  augmented  measurement  vector  defined  as: 

X 
Zn 

•  is  the  (5  X  4)  augmented  regressor  defined  as: 

..  f  / 1 


•  K  is  the  (5  X  1)  augmented  measurement  noise  vector  defined  as: 


W 

Wn 
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The  goal  is  to  derive  a  Kalman-like  update  equation  to  refine  the  unknown  GPS 
parameters  vector  estimate,  {X),  and  its  covariance  matrix,  (Px),  both  produced  by 
the  closed-form  algorithm  developed  in  Section  3.1.  In  the  derivation  that  follows, 
to  distinguish  the  estimates  X  and  Px  as  produced  by  the  closed-form  algorithm 
from  the  estimates  that  will  be  obtained  through  the  Kalman  update,  the  following 
notation  is  used: 

•  X  and  Px  represent  the  estimate  and  the  estimation  error  covariance  prior 
to  the  update;  and 

•  X  and  Pf^  represent  the  estimate  and  the  estimation  error  covariance  follow¬ 
ing  the  update. 

In  order  to  obtain  the  updated  estimates  from  the  augmented  linear  regression  in 
Equation  (3.40),  it  is  necessary  to  derive  the  covariance  of  the  augmented  noise  vector 
Va-  Since  the  statistics  of  the  noise  components  in  14  have  already  been  determined, 
the  equation  error  covariance  matrix,  Ra,  is  derived  as  follows: 


R 


a 


=  E{} 

Iv,^} 

=  E{ 

W 

[W^  Wn]} 

.  . 

WW’^  WWn 

=  E{ 

WnW^  wl 

Px  P 


(3.41) 


The  updated  minimum  variance  GPS  solution  estimate  and  associated  covari¬ 
ance  matrix  are  given  by  the  expressions: 


k*  =  (3.42) 

Pt  =  (3.43) 
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The  expressions  in  Equation  (3.42)  and  Equation  (3.43)  are  sufficient  to  obtain  the 
required  updates  but  it  is  desirable  to  manipulate  and  reduce  the  equations  into 
the  more  familiar  and  computationally  efficient  form  of  the  classical  Kalman  filter 
update  equations.  The  expressions  can  be  rewritten  in  expanded  form  as; 


[I  h] 

[[/  h] 


1 _ 

P 

-1 

X 

(t2 

Zn 

- 1 

\ 

P 

1 

I 

T 

P 

0-2 

-1 


(3.44) 


(3.45) 


To  assist  in  the  reduction  of  the  expressions,  the  parameters  are  redefined  by 
scaling  them  as  follows; 


=!■ 

P*  ^  a-^Pi 

p  ^  cr^p 


In  terms  of  the  redefined  parameters,  the  expressions  become; 


'x  =  Pi\ih\ 

1 

1 _ 

-1 

—*■ 

X 

p^  1 

Zn 

P^ 


[[/  h] 


'  Px 

P 

-1 

I 

p^ 

1 

hT 

The  inverse  of  the  redefined  noise  covariance  matrix  is  calculated  by  defining 
a  piecewise  matrix  inverse  and  algebraically  determining  the  inverse  of  each  piece; 

Px  P 

p^  1 


—  X 

A 

a 

— 

T 

a 

c 

(3.46) 
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From  Equation  (3.46),  it  is  obvious  that, 


1 _ 

A  a 

I  0 

1 - 

1 _ 

T 

a  c 

0  1 

(3.47) 


The  following  four  equations  are  obtained  by  piecewise  equating  the  members  of  the 
matrix  on  the  left  hand  side  of  Equation  (3.47)  to  those  on  the  right  hand  side: 


•  A  +  paF  =  / 

•  A  +  =  0 

•  Pxo,  +  cp  =  0 

•  p^a  +  c  =  1 

From  the  four  equations  above,  expressions  for  A,  a,  and  c  are  determined  alge¬ 
braically  in  terms  of  only  Px  and  p. 


A  =  {Px-pp'^r^ 

a  =  -{Px-PP^y^P 
c  =  1 + p'^iPx  -  pp^y^p 


(3.48) 

(3.49) 

(3.50) 


With  A,  a,  and  c  known,  it  is  possible  to  describe  ^ .  Further  to  the  above  equations 
the  following  relationships  are  noted: 


a  =  -Ap  (3.51) 

c  =  1  +  p^  Ap  (3.52) 


Using  Equations  (3.51)  and  (3.52),  the  following  two  expressions  that  will  be  ex¬ 
tremely  useful  in  the  simplification  of  Equations  (3.44)  and  (3.45)  are  formed: 

A  ha^  =  A  —  hp^A  =  {I  —  hp^)A 
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a-+  ch  =  —Ap  +  h  +  hp^Ap  =  h  —  [I  —  hp'^)Ap 


The  expressions  in  Equations  (3.44)  and  (3.45)  can  be  rewritten  by  substituting 
R~^  by  its  equivalent  form  as  defined  in  Equation  (3.46)  and  simplified  by  multiplying 
out  the  matrices. 


1++  1 

A  a 

X 

X  =  Pill  A] 

T 

a  c 

Zn 

=  [(A  T  ha^^X  -(-  (fl  ch^Z^i] 

=  Px{A  +  ha'^)X  +  Px{a  +  ch)Zn 

=  P^{I-hp'^)Ai  +  P^[h  -  {I  -  hp'^)Xp\Zn  (3.53) 


A  a 

1 

T 

a  c 

I 

1 _ 

=  [A  +  chh^  +  ha^  +  ah^]~^ 

=  [A  +  (a  +  ch)h^  +  ha^]~^ 

=  [A  +  {h  —  {I  —  hp^)Ap)h^  +  ha^]~^ 

=  [A  +  hh^  ~  (1  “  hp^)Aph^  +  ha^]~^ 

=  [A  +  hh^  —  Aph^  +  hp^Aph^  —  hp^A]~^ 

=  [hh^  +  ~  hp^)A  —  {I  —  hp^)AphF]~^ 

=  [hh^  +  (J  -  hp^)A{I  -  ph^)]~^  (3.54) 


The  inverse  of  the  (I  —  ph^)  matrix  in  the  P^  expression  in  Equation  (3.54)  is 
calculated  by  applying  the  Matrix  Inversion  Lemma  (MIL): 

(/  -  =  I +  h{l-  p^h)~^p^  =  /  +  ^ 
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By  using  the  above  result  and  introducing  the  expression  for  A  in  Equation  (3.48), 
Equation  (3.54)  becomes: 


1  1  ph^)-T' 


P+  =  [hh  +il  +  ^^hp)-  A{I+^^ 


=  +  (J  4:.  ^  fep^)  \Px-PP^)  ^(J  +  ^  )  1] 


To  simplify  the  expression  for  further,  it  is  necessary  to  define  a  new  inter¬ 
mediate  variable,  Y ; 


(3.55) 


The  expression  for  Y  defined  in  Equation  (3.55)  is  further  simplified  as  follows: 


Y  =  (/  + 


1  ,  'Tv  /  y  .  1 


1  —  p'^h 
1 


hp  ){Px  -PP  + 


Pxhp'^- 
1 


1  —  p^h 


1  —  p^h 
-uJT 


p^hpp^) 


pTh^ 


p^h 


hPPyh  —  hPp  +  p^h  —  1  j. 

=  Pii  +  — ^  Iri;, - PP  + 


1 


(1  —  p'^hy 


1  —  p'^h 


(Pxhp^  +  ph'^  Px) 


=  Px  + 


hFPxh-1 


PP  + 


1 


{1-pThy^^  ‘  i-pTh 


■{Pxhp^  +  ph'^  Px) 


(3.56) 


By  introducing  the  expression  for  Y  defined  in  Equation  (3.55)  and  applying 
the  MIL,  the  expression  for  P^  is  further  simplified: 


P^ 


=  (Y-^  +  hh^)-'^ 

=  Y  -Yh{l  +  h'^Yh)-'^h'^Y 

1 

=  Y _ _ _ Ykh^Y 

lAhTYh 


(3.57) 
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It  is  noted  from  the  simplified  expression  for  Px  given  in  Equation  (3.57)  that  no 
matrix  inversion  operation  is  required  for  the  update  of  the  covariance  matrix  since 
both  Equations  (3.56)  and  (3.57)  do  not  involve  any  matrix  inversion. 

Attention  must  now  be  turned  towards  the  update  equation  for  the  GPS  solu- 

A 

tion  estimate  (X).  In  Equation  (3.53)  the  term,  {P^{I  —  hp^)A},  is  algebraically 
manipulated  to  allow  for  further  simplification  of  the  equation.  The  manipulation  in¬ 
volves  using  the  result  obtained  in  the  covariance  update  equation  in  Equation  (3.57): 


PUI  -  hi')A 


{I  -  P^hh^){I  -  ph^)-^ 


_  h^p 
1  —  p^h 

P^hh^ 


1  + 


1 


1  —  p'^h 


{p  ~  Px^)^ 


P^hh^ 


(3.58) 


The  expressions  in  Equation  (3.58)  and  Equation  (3.48)  are  replaced  into  Equa¬ 
tion  (3.53)  and  the  equation  is  further  manipulated  to  obtain: 

=  P^{I  -  hp'^)Ak~  +  P^[h  -  {I  -  hp^)Ap]Zn 
=  P^{I  -  hp^)Ak~  +  P^hZn  -  P^I  -  hp^)ApZn 
=  p+(/  -  hp^)Ak~  +  p+hZr,  -  (/  -f  y:^{p  -  p^h)h'^)pZn 

=  P+(7  -  hp'^)Ak~  +  P+hZn  -  pZn  -  Y^^{p  -  P^h)h^pZn 

=  P+{I-hp^)Ak~  +  P+hZn-pZn-:^^^^{p-P^h)Z^ 

=  P+{I  -  hp^)Ak~  +  P^hZ^  -  pZn  -  Y^^ipZn)  +  J^^iP^hZr,) 

=  p${i  -  hp^)Ak~ + p^hz^  - 
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"  I- ft  -  J  _  pTftO’^J  +  1  _  pT^(^x'‘2«) 

=  +  riW'"  - 

=  ^'  +  riW'’^'*’’^'  ■ 

=  i'  +  Y3^(Pl(.-p)(Z„-rf“)  (3.59) 

By  comparing  the  expression  in  Equation  (3.59)  to  the  classical  Kalman  filter 
estimate  update  equation,  it  is  noted  that  by  defining  a  filter  gain  K: 

Equation  (3.59)  is  reduced  to  the  form  of  the  classical  estimate  update  equation  and 
it  can  be  rewritten  as: 

=  JC~  +K{Zn-h^X  )  (3.60) 

Although  Equation  (3.60)  has  the  same  form  as  that  of  the  classical  Kalman 
filter,  it  must  be  noted  that  the  filter  gains  K  are  not  the  same.  The  main  difference 
is  due  to  the  introduction  of  the  correlation  p  between  the  new  measurement  and 
the  previous  estimate  which  in  the  case  of  the  classical  Kalman  filter  is  zero.  As 
is  the  case  for  the  Kalman  filter  update  "equations,  an  equivalent  expression  for  the 
filter  gain  K  must  be  obtained  in  terms  of  the  covariance  prior  to  the  update  (Px)- 
Using  the  previous  result  obtained  for  in  Equation  (3.57),  it  is  most  convenient 
to  express  the  filter  gain  K  in  terms  of  Y  as  defined  in  Equation  (3.55),  which  is  a 
function  of  {Px)- 
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1 


1 


1  —  p^h 
1 


[Yh- 


1  +  K^Yh 
1 


1  —  p^h 

— I — r - - - 

1-pThh  +  hTYh 


1  +  h^Yh 

Yh-p] 


Yhh^Y)h  -  p] 
Yh{h^Yh)  -  p] 


(3.61) 


Using  the  new  expression  for  K  obtained  in  Equation  (3.61),  the  covariance 
update  equation  (  Equation  (3.57)  )  must  be  rewritten  in  terms  of  the  filter  gain  K. 
By  manipulating  Equation  (3.61),  the  following  expression  is  obtained: 


(1  -p^h)K  +  p  = 


1 

1  +  h^Yh 


Yh 


This  expression  is  inserted  into  Equation  (3.57)  to  obtain  the  covariance  update 
equation  in  the  desired  form: 


P+  =  Y  -[{l-p'^h)K  +  p]h^Y 

=  {I-  [(1  - p^h)K  +  p\h'^]Y  (3.62) 

The  derivation  of  the  Kalman-like  update  equations  in  the  desired  form  is  now  com¬ 
plete. 

To  summarize  the  results  obtained,  the  final  Kalman  update  equations  as  pre¬ 
sented  in  Equation  (3.56),  Equation  (3.61),  Equation  (3.60),  and  Equation  (3.62) 
are: 


Y 

K 

X 


Px  + 


hFPxh-l 


(1  —  p^hy 

^  1  ^ 


l-p^hn  +  K^Yh 


pp  + 


Yh 


1 


{Pxhp^  P pK^' Px) 


To-y 


1  —  p^h 

■v\ 


=  k  +K{Zn-hyk  ) 

=  {I-[{l-p^h)K  +  p]h^}Y 
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The  above  equations  are  used  in  the  Matlab  implementation  of  the  Kalman  up¬ 
date  experimental  algorithm.  The  Kalman  update  algorithm  is  intended  to  refine  the 
GPS  solution  estimate  in  a  direct  and  non-iterative  manner.  Although  this  can  be 
achieved  in  most  cases,  it  was  determined  experimentally  that  sometimes  a  second 
application  of  the  algorithm  is  required  to  obtain  a  better  solution  estimate.  Re¬ 
calling  that  the  new  measurement  used  by  the  Kalman  update  algorithm  is  actually 
the  pseudorange  equation  in  Equation  (3.3)  that  has  been  linearized  about  the 
position  estimate  produced  by  the  closed-form  algorithm,  implies  that  how  well  the 
linearization  fits  the  true  unknown  GPS  parameters  is  dependent  on  how  good  the 
solution  produced  by  the  closed- form  algorithm  is  to  start.  In  order  to  alleviate  this 
undesired  dependency,  after  the  Kalman  Update  algorithm  has  been  applied  once, 
and  produces  an  improved  solution  estimate.  Equation  (3.3)  is  once  again  linearized 
about  the  improved  position  estimate  producing  a  better  new  measurement.  The 
Kalman  update  algorithm  is  applied  a  second  time  keeping  in  mind  that  the  estimate 
prior  to  the  update  {X  )  and  the  covariance  prior  to  the  update  {Px)  are  the  esti¬ 
mates  produced  by  the  closed-form  algorithm,  not  the  solution  obtained  as  a  result 
of  the  previous  application  of  the  Kalman  update.  Theoretically,  this  process  can 
be  continued  recursively  until  convergence  to  the  best  possible  solution  is  achieved; 
however,  it  was  shown  experimentally  that  after  the  second  application  of  the  algo¬ 
rithm,  the  change  in  the  solution  estimate  is  insignificant.  Consequently,  continued 
iterations  are  not  required;  the  algorithm  needs  to  be  applied  just  twice.  The  first 
application  is  strictly  to  obtain  a  suitable  position  estimate  about  which  to  perform 
a  valid  linearization  of  the  pseudorange  equation  and  the  second  application  is 
to  calculate  the  final  GPS  solution  estimate  and  its  covariance. 

A  quick  verification  of  the  derived  update  equations  in  Equations  (3.56),  (3.61), 
(3.60),  and  (3.62)  is  carried  out  to  confirm  the  validity  of  the  new  equations.  The 
special  case  of  the  classical  Kalman  filter  where  p  =  0  is  considered.  For  this  spe¬ 
cial  case  of  no  correlation,  the  well  known  classical  Kalman  update  formulae,  are 
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recovered; 


V 


K 


X 


Px 

k~  +K{Zn-h'^k~) 
(/  -  Kh^)Px 


3.4  Alternate  Stochastic  Closed-form  Solution 

This  section  presents  the  mathematical  derivation  of  an  alternate  stochastic 
closed-form  solution  that  is  closely  related  to  the  closed-form  algorithm  developed 
in  Section  3.1.  This  alternate  closed-form  algorithm  attempts  to  alleviate  some  of 
the  undesirable  characteristic  associated  with  the  closed-form  algorithm  developed  in 
Section  3.1.  The  development  of  this  algorithm  is  presented  strictly  to  demonstrate 
a  potentially  viable  approach  in  producing  an  alternate  stochastic  closed-form  solu¬ 
tion.  The  main  difference  between  this  algorithm  and  the  previous  one  is  that  the 
differencing  of  the  n  expanded  pseudorange  equations  presented  at  Equation  (3.1) 
will  be  done  in  a  cyclical  fashion.  The  concept  of  cyclical  differencing  is  based  on 
the  difference  linearization  approach  demonstrated  by  Krause  [14]  for  the  develop¬ 
ment  of  a  deterministic  closed-form  solution.  The  solution  in  [14]  uses  differencing 
to  formulate  a  computation  basis  from  which  the  user  position  is  backed  out  and 
uses  a  nonlinear  auxiliary  equation  to  compute  the  user  clock  bias;  whereas,  the  goal 
in  this  development  is  to  form  a  linear  regression  consisting  of  n  equations  in  four 
unknowns,  three  user  position  coordinates  and  the  user  clock  bias. 

The  cyclical  differencing  scheme  involves  subtracting  the  (i-|-l)‘^  equation  from 
the  equation  for  i  =  1, 2, ...,  n  —  1  where  n  is  the  number  of  available  pseudorange 
measurements.  The  last,  (n*^),  linearized  equation  is  formed  by  subtracting  the  first 
equation  from  the  n**  equation.  The  most  noteworthy  motivation  to  attempt  this 
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development  is  that  the  cyclical  differencing  will  eliminate  the  nonlinear  terms  while 
retaining  n  linearized  equations,  unlike  the  previous  closed-form  solution  presented 
in  Section  3.1  which  produced  (n  —  1)  linearized  equations  and  was  left  with  the 
nonlinear  pseudorange  equation  for  use  as  an  auxiliary  equation.  It  was  initially 
anticipated  that  this  alternate  closed-form  solution  would  be  capable  of  producing 
estimates  of  the  four  GPS  solution  parameters  with  only  four  pseudorange  measure¬ 
ments,  unlike  the  closed-form  algorithm  presented  in  Section  3.1  which  requires  a 
minimum  of  five. 

3.4’ 1  New  Linear  Regression  in  n  Equations.  The  initial  portion  of  the 
derivation  of  this  alternate  closed-form  solution  to  the  pseudorange  equations  is  iden¬ 
tical  to  that  presented  in  Section  3.1.  A  natural  starting  point  for  this  development 
is  the  point  at  which  the  developments  differ.  Consequently,  the  development  of  the 
alternate  closed-form  solution  begins  with  Equation  (3.1); 

-  6^  -  2xiUx  -  2yiUy  -  2ziUz  -t-  2Rib  = 

R]  -  xf  -  y1  -  zf  -  2RiWi  +  2bwi  -|-  w1 

The  cyclical  differencing  is  applied  resulting  in  n  equations  that  are  linear  in  the 
unknown  variables  and  can  be  expressed  as: 

xffUx  T  (2/2+1  T  (-^i+l  T  {.^Ri  Ri.\-l^b 

^{Ri-Rf+i+x^i+i-x^i+yi+i-yf+z^+^-zf)-\-Ri+iWi+x-RiWi-^bwi-bwi+x-\-^{w^-w^^.^) 

(3.63) 

Equation  (3.63)  is  understood  to  encompass  n  equations  for  i  =  l...n,  where  it  is 
understood  that  the  first  equation  also  serves  as  the  (n  -|- 1)*^  equation. 

To  simplify  Equation  (3.63),  the  approximation  at  Equation  (3.5)  is  introduced 
as  was  done  in  the  previous  derivation  and  Equation  (3.63)  becomes: 
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,(^z+l  “i"  (^z+l  Ui^'^y  ('^2+1  “1"  -^2+1  )^  ^ 

^{Ri-R'f+i+x^i+i-x^i+yi+i-yH^hi-^i)+i^-^)i^i-^i+^)+\(^^-'^^+i)  (3-^4) 

By  redefining  the  following  terms  : 

•  (6  -  R){wi  -  Wi+i)  +  l{wf  -  =»  K-; 

•  K^i  -  ^1+1  +  ^i+i  -x]  +  yf+i  -  yf  +  ^-+1  -  zf)  ^  Zi 

Equation  (3.64)  reduces  to  a  linear  regression  form  of  n  equations  in  four  unknowns 
expressed  as: 


{xi^i  -  Xi)ux  +  (i/i+i  -  yi)uy  +  {zi+i  -  Zi)uz  +  {Ri  -  Ri^\)h  ft!  Zi  +  Vi  (3.65) 


At  this  point,  the  linear  regression  in  Equation  (3.65)  can  be  written  in  matrix 
notation  form  as: 

Z  =  HX  +  V,  (3.66) 

where 


Z  is  the  measurement  vector  of  dimension  n  given  by: 
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•  H  is  the  (n  x  4)  regressor  matrix  given  by: 

X2  -Xl  J/2  -  yi  Z2  -  Zi  Ri- 

X3  -X2  t/3  -  2/2  ^3  -  Z2  R2  —  Rz 

H=  .  .  .  , 

Xn  Xfi—\  yn  yn — 1  Z^i  R-fi—X  Rn 

Xl  Xfi  y\  2/n  ■^1  Rn  Rl 

•  X  is  the  vector  of  unknowns  given  by: 

y'x 

Uy 

Uz 

b 

•  y  is  the  equation  error  vector  of  dimension  n  given  by: 

'  y 
y 

y=  : 

_K 

It  must  be  noted  that  the  the  (n  x  4)  regressor  matrix  (H)  is  rank  deficient,  viz., 
rank{H)  =  3;  therefore,  this  approach  is  incapable  of  producing  a  GPS  solution  in 
the  deterministic  case. 

3.4.2  Noise  Statistics.  Using  the  linear  regression  obtained  in  Equa¬ 
tion  (3.66),  it  is  possible  to  obtain  an  estimate  of  X  if  the  statistics  of  the  equation 
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error  vector  V  are  known.  The  statistics  of  V  must  be  derived  from  the  known 

— f 

statistics  of  the  pseudorange  measurements  noise,  Wi.  The  equation  error  vector  V 
is  known  to  be  zero  mean;  therefore,  only  the  covariance  matrix  of  the  equation 
error  vector  V  needs  to  be  derived.  The  covariance  matrix  is  an  (n  x  n)  symmetrical 
matrix  whose  entries,  in  the  most  general  form,  are  given  by: 


E{ViVj}  =  E{[(b-R)(wi  -  Wi+i)  +  ]^{wj  -  u;-+i)][(6-  R){wj  -  Wj+t)  +  -  w]^^)]} 

(3.67) 

The  subscripts  i  and  j  in  Equation  (3.67)  take  on  values  from  1  through  n  and 
it  is  understood  that  a  subscript  of  (n  +  1)  is  equivalent  to  a  subscript  of  1.  There 
are  three  separate  cases  of  covariance  entries  that  must  be  considered: 


1.  i  =  j.  These  entries  represent  the  diagonal  elements  of  the  covariance  matrix 
and  are  annotated  as  Pvd- 

2.  i  =  j  +  I  OT  j  =  z  +  l.  These  entries  represent  the  adjacent  elements  and 
populate  the  matrix  cells  adjacent  to  the  diagonal  elements.  Matrix  elements 
(l,n)  and  (n,  1)  are  also  adjacent  elements  since  it  is  understood  that  (n  +  1) 
is  equivalent  to  1.  These  adjacent  elements  are  annotated  as  Py^. 

3.  i  j  1  and  j  ^  i  +  1  and  i  ^  j.  These  entries  represent  all  the  remaining 
elements  of  the  covariance  matrix  and  are  annotated  as 

The  noise  covariance  matrix  Py  whose  elements  are  described  by  E{ViVj}  is  a 
symmetrical  (n  x  n)  matrix  expressed  as: 


Pvd 

Pva 

Pvr 

Pv. 

Pvd 

Pva 

Pvr 

**  • 

Pvr 

Pvr 

.  Pva 

Pvr 

Pvr 

Pva 

Pvr 

■■  ^ 

■ . 

Pvr 

Pva 

Pvd 

Pva 

Pvr 

Pva 

1 _ 

(3.68) 


3-38 


By  recognizing  that  the  expectation  operator  (E)  is  a  linear  operator,  the  ex¬ 
pression  for  E{ViVj}  in  Equation  (3.67)  can  be  expanded  as  follows; 

E{ViVj}  =  E{{b- Rf  {wi-wi+i){wj -Wj+i)}  +  E{^{b- R) 

[(u)i  -  Wi+i){w]  -  +  {wj  - 

+E{^{w^  - 

=  (6  -  RfEKwiWj  -  WiWj+i  -  Wi+iWj  +  tOj+iUJj+i)}  -f  -  R) 

E{{Wiw]  -  -  Wi+iw]  -1-  Wi+iw]^^  +  WjWi  -  Wjwf^^  -  Wj+iwf 

+u;j+i«)?+i)}  +  ^E{{w^w]  -  (3.69) 

From  Equation  (3.69)  it  is  possible  to  find  the  equation  error  statistics  by  using  the 
defined  elemental  noise  statistics  of  the  pseudorange  measurement  noise,  (tn).  In 
order  to  proceed,  each  of  the  three  cases  of  the  covariance  matrix  entries  must  be 
considered  individually,  where  Equation  (3.69)  is  further  reduced  for  each  case. 

3. 4- 2.1  Noise  Covariance,  Diagonal  Elements.  The  diagonal  ele¬ 
ments  of  the  covariance  matrix  are  represented  by  the  case  where  j  =  i.  Making  this 
substitution  into  Equation  (3.69)  reduces  the  expression  for  the  diagonal  covariance 
elements  to: 

Pvd  =  (b-  R)'^E{{wf  -  2wiWi+i  +  ro'+i)} 

+{b  -  R)E{{w^  -  wiwf^^  -  Wi+iWi  +  u;-+i)} 

-  ^t+i)} 

=  {b-Ry{(T^-0  +  (T^)  +  {b-R){0)  +  ^i3a^-2a^  +  Za^) 

=  2(6  —  RYcr^  -t- 

3. 4- 2.2  Noise  Covariance,  Adjacent  Elements.  The  adjacent  elements 
of  the  covariance  matrix  are  represented  by  the  case  where  j  -|- 1  =  i.  Making  this 


3-39 


substitution  into  Equation  (3.69)  reduces  the  expression  for  the  adjacent  covariance 
elements  to: 


Pva  -  (b- RyE{{wiWj -wf  -  Wi+iWj +Wi+iWi)} +  ^{b- R) 

E{{wiw]  -  twf  -  Wi+iw]  +  Wi+iW^  +  Wjwf  - 
+\E{{wfw]  -  wf  -  +  wf+^wf)} 

=  (6  —  ^)^(0  —  0  —  cr^  +  0)  +  ^(6  —  -R)(0)  +  —  3cr^  —  0-^  +  cr^) 

=  _(6-B)V^_icr‘  = 

3.4-2. 3  Noise  Covariance,  All  Remaining  Elements.  Equation  (3.69), 
in  its  current  form,  represents  the  case  for  all  the  remaining  covariance  matrix. 

Pvn  =  (6  -  .R)^(O)  +  i(6  -  .R)(0)  +  i((T^  -  -  <7^  +  <T^)  =  0 

Since  the  individual  covariance  matrix  elements  Pvn,  Pvai  and  Pvr  have  been 

— ♦ 

determined,  the  covariance  of  the  equation  error  vector  V  can  be  expressed  as: 

Py^a\^a^A{b-Rf)Rn,  (3.70) 

where  ^  is  a  (n  x  n)  matrix  whose  diagonal  elements  are  2,  the  elements  adjacent 
to  the  diagonal  elements  are  —1,  and  all  the  remaining  elements  are  zero.  As  an 
example,  the  Rn  matrix  for  n  =  5  is  a  (5  x  5)  matrix  defined  as  follows: 

2-10  0-1 
-12-10  0 
Rs=  0-12-10 
0  0-12-1 
-10  0-12 
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An  interesting  characteristic  of  the  Rn  matrix  is  that  it  is  always  positive 
semi-definite.  Rn  always  has  exactly  one  eigenvalue  equal  to  zero;  hence,  it  is  not 
invertible.  Recalling  the  solution  for  the  minimum  variance  estimate  of  the  unknown 
parameter  vector  X  derived  in  Equation  (3.16),  the  inverse  of  is  used  as  a  weight¬ 
ing  matrix.  Since  is  not  invertible,  the  unknown  parameter  estimation  is  obtained 
using  the  following  equation  where  the  inverse  of  is  substituted  by  the  generalized 
inverse  of  Rn‘. 

k  =  (H'^rIu)-^IFrIz,  (3.71) 

where  ill  is  the  generalized  inverse  of  Rn- 

To  get  a  deeper  appreciation  of  the  problem  at  hand,  which  is  indeed  more 
complex  than  it  appears,  the  concept  of  the  generalized  inverse  is  examined  closely. 
There  are  a  number  of  possible  approaches  to  obtaining  a  generalized  inverse  of  a 
matrix,  but  given  the  symmetric  characteristic  of  the  matrix,  using  the  Singular 
Value  Decompostion  (SVD)  approach  is  the  most  convenient  [4].  Performing  a  SVD 
on  a  real  symmetric  produces  the  following  result: 

Rn  =  USU^  (3.72) 

where  5  is  a  (n  x  n)  diagonal  matrix  whose  diagonal  elements  are  simply  the  singular 
values  or  eigenvalues  of  Rn  and  17  is  a  (n  x  n)  unitary  matrix.  Using  the  SVD  results, 
the  generalized  inverse  oi  Rn  is  simply  given  by: 

rI  =  (3.73) 
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where  the  generalized  inverse  of  S  is  given  by: 

t  0  .  0 

0  ±  0  ...  0 

0  •••  0 

0  •••  0  0 

0  .  0 

The  eigenvalues  of  Rn  are  represented  by  Aj  but  it  must  be  reemphasized  that  Rn 
has  one  eigenvalue  (A^  =  0). 

The  rank  deficiency  is  now  evident;  furthermore,  the  last  column  of  the  unitary 
matrix,  U,  is  arbitrary  and  insignificant.  The  impact  of  this  rank  deficiency  is 
that  although  the  linear  regression  in  Equation  (3.66)  appears  to  be  of  order  n, 
it  is  actually  of  order  (n  —  1).  As  was  the  case  in  the  closed-form  algorithm,  a 
minimum  of  five  pseudorange  measurements  are  required  to  produce  a  solution  using 
this  alternate  closed-form  approach.  To  avoid  the  apparently  excessive  order  of  the 
linear  regression  and  the  requirement  to  work  with  a  generalized  inverse,  the  linear 
regression  in  Equation  (3.66)  can  be  transformed  into  an  equivalent  form  of  order 
(n-1). 

A  transformation  matrix,  U^,  is  defined  and  recalling  that  the  last  column  of 
U  is  insignificant  allows  for  the  transformation  matrix  to  be  partitioned  as  follows: 

T 

U^=  .. 

t 

where  t  is  the  last  row  of  the  matrix.  A  new  linear  regression  is  obtained  by 
premultiplying  the  linear  regression  in  Equation  (3.66)  by  the  ((n  —  1)  x  n)  trans¬ 
formation  matrix,  T : 
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TZ  =  TEX  +  TV,  (3.74) 

A  simple  linear  regression  is  formulated  by  redefining  the  terms  in  Equation  (3.74) 
as  follows: 

•  TZ  ^  Zt 

•  TH^Ht 

.  ry  ^ 

The  transformed  linear  regression  is  represented  by  : 

Zt  =  HtX  +  Vt,  (3.75) 


where 

•  is  the  transformed  measurement  vector  of  dimension  (n  —  1) 

•  Ht  is  the  transformed  regressor  matrix  of  dimension  ((n  —  1)  x  4) 

•  Vt  is  the  transformed  equation  error  vector  of  dimension  (n  —  1) 

The  covariance  matrix  of  the  equation  error  vector  Vy  is  now  calculated  as  follows: 


E{VtV^}  =  E{TW(TWf} 
=  E{TWW^T'^} 
=  TE{WW^}T'^ 


=  TRnT^ 


(3.76) 


The  expression  for  the  noise  covariance  obtained  in  Equation  (3.76)  is  recog¬ 
nized  as  the  diagonal  S  matrix  obtained  from  the  SVD  operation  in  Equation  (3.72) 
without  its  ri^h  row  and  n*h  column  of  zeros.  The  covariance  matrix  can  be  defined 
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as  E{VtV^]  =  Sn-\,  where  5„_i  is  given  by: 


1 


Ai  0  •  •  •  0 

0  A2  0  : 

:  •••  0 

0  •  •  •  0  A^-i 


The  diagonal  form  of  the  transformed  noise  covariance  simplifies  the  task  of  calcu¬ 
lating  its  inverse  which  is  explicitly  given  by: 


S. 


-1 

n—l 


i  0  •••  0 

0*0  = 

;  •••  0 


The  minimum  variance  estimate  of  the  unknown  GPS  parameter  vector  based  on 
the  transformed  linear  regression  in  Equation  (3.75)  can  thereby  be  calculated  using 
the  following  expression: 


X  =  (3.77) 

Equation  (3.77)  is  the  most  convenient  for  use  in  the  implementation  of  the 
experimental  algorithm  in  Matlab.  Preliminary  experimental  trials  on  the  alternate 
closed-form  algorithm  demonstrated  that  the  use  of  a  supplementary  algorithm  is 
required  to  refine  the  solution  produced  by  the  alternate  closed-form  algorithm,  as 
was  the  case  for  the  closed-form  algorithm.  Both  the  maximum  likelihood  and  the 
Kalman  update  algorithms  can  be  used  to  supplement;  however,  now  there  is  no 
obvious  choice  of  auxiliary  equation  to  use.  Any  of  the  original  n  pseudorange 
equations  can  be  used  as  the  auxiliary  equation.  The  application  of  the  maximum 
likelihood  to  the  alternate  closed-form  solution  is  identical  to  that  which  was  used  for 
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the  closed-form  solution.  In  order  to  apply  the  Kalman  update  algorithm  it  would 
be  necessary  to  determine  the  correlation  between  the  solution  produced  by  the 
alternate  closed-form  algorithm  and  the  selected  auxiliary  equation.  This  correlation 
is  not  the  same  as  the  one  determined  for  the  closed-form  solution;  furthermore,  the 
development  is  more  involved  due  to  the  more  complex  form  of  the  noise  covariance 
matrix.  The  development  of  the  Kalman  update  algorithm  tailored  for  application 
to  the  alternate  closed-form  solution  is  not  addressed  in  this  thesis. 

3.5  Conclusion  of  Mathematical  Derivations 

This  chapter  presented  the  derivations  of  the  four  algorithms  that  are  used  in 
the  Matlab  implementations  for  experimental  analysis.  The  four  algorithms  are; 

1.  The  Closed-form  algorithm 

2.  The  Maximum  likelihood  algorithm 

3.  The  Kalman  update  algorithm 

4.  The  Alternate  Closed-form  algorithm 

The  first  three  algorithms  axe  thoroughly  analyzed  experimentally  in  Chap¬ 
ter  4;  however,  only  the  Closed-form  algorithm  developed  in  Section  3.1  provides  a 
stand-alone  GPS  solution  estimate,  the  two  latter  are  used  to  refine  the  solution  pro¬ 
vided  by  the  closed-form  algorithm.  The  maximum  likelihood  algorithm  developed  in 
Section  3.2  accomplishes  the  task  of  refining  the  solution  but  it  inherently  possesses 
two  serious  drawbacks: 

•  It  depends  on  optimization  techniques  that  are  recursive  in  nature  hence  the 
closed-form  appeal  is  lost. 

•  It  does  not  have  the  capability  to  provide  the  estimation  error  covariance  read¬ 
ily. 
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As  a  result  of  these  drawbacks,  the  more  elegant  Kalman  update  supplementary 
algorithm  was  developed. 

The  Kalman  update  algorithm  maintains  a  closed-form  approach  despite  the 
fact  that  a  portion  of  the  algorithm  must  be  applied  twice  to  ensure  a  good  solution. 
The  algorithm  also  retains  the  capability  to  calculate  the  estimation  error  covariance. 
Based  strictly  on  the  mathematical  attributes  of  the  Kalman  update  algorithm,  it  is 
the  preferred  supplementary  algorithm.  The  performance  of  each  of  the  algorithms  is 
assessed  experimentally  to  provide  a  broader  comparison  based  also  on  performance. 

The  alternate  closed-form  algorithm  based  on  cyclical  differencing  linearization, 
was  expected  to  provide  a  more  mathematically  aesthetic  and  balanced  solution  than 
the  closed-form  solution.  The  algorithm  produces  a  singular  regressor  matrix,  H, 
preventing  this  approach  from  being  used  in  the  deterministic  case.  As  well,  in 
the  stochastic  case,  the  algorithm  produces  a  positive  semi-definite  noise  covariance 
matrix  resulting  in  an  undesirable  rank  deficiency;  consequently,  just  like  the  closed- 
form  algorithm,  it  still  requires  at  least  five  pseudorange  measurements  in  order  to 
produce  a  GPS  solution.  The  performance  of  the  alternate  closed-form  algorithm 
needs  to  be  compared  to  that  of  closed-form  algorithm  through  experimental  analysis 
to  confirm  which  closed-form  solution,  if  any,  is  preferable. 

The  most  interesting  realization  that  was  brought  to  light  through  the  mathe¬ 
matical  development  of  the  two  stochastic  closed-form  solutions  is  the  set  of  conse¬ 
quences  of  working  with  stochastic  pseudorange  equations.  In  both  solutions,  the  av¬ 
erage  pseudorange  approximation  at  Equation  (3.5)  is  required  in  order  to  obtain  an 
expression  for  the  noise  covariance  in  which  the  dependency  on  any  of  the  estimation 
parameters  can  be  eliminated.  This  approximation  is  not  required  in  a  deterministic 
development.  It  turned  out  that  both  solutions  required  at  least  five  pseudorange 
measurements  to  produce  a  GPS  solution.  In  the  closed-form  algorithm  this  was 
a  direct  result  of  the  methodology  used  in  differencing  the  pseudorange  equations; 
however,  the  alternate  closed-form  algorithm  was  expected  to  resolve  this  problem 
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but  was  unable  to  do  so  due  to  the  positive  semi-definite  noise  covariance  matrix, 
with  a  single  zero  eigenvalue  guaranteed,  obtained.  In  conclusion,  when  accounting 
for  measurement  noise,  it  is  impossible  to  obtain  a  closed-form  fix  using  linear  math¬ 
ematics  only;  furthermore,  a  direct  (closed-form)  solution  can’t  be  obtained  using 
only  four  pseudorange  measurements. 
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IV.  Experimental  Results  and  Analysis 

This  chapter  presents  the  experimental  portion  of  the  thesis.  The  first  part 
discusses  how  the  experiment  was  set  up  and  how  the  Monte  Carlo  trials  were  run. 
Next  is  presented  a  brief  discussion  of  how  the  results  are  intended  to  be  interpreted. 
The  experimental  results  are  then  presented  and  the  chapter  sums  up  with  a  detailed 
analysis  of  the  results. 

4.1  Overview 

The  experimental  results  encompass  the  outcome  from  two  test  environments 
which  were  used  to  conduct  the  experimental  trials.  The  first  test  environment 
consists  of  twelve  simulated  NAVSTAR  GPS  satellite  scenarios  chosen  to  represent 
conditions  typically  encountered  in  the  real  world.  The  twelve  satellite  scenarios  were 
generated  using  the  GPSoft  Satellite  Navigation  Toolbox  for  Matlab  [10].  The  second 
test  environment  consists  of  three  scenarios  of  pseudo-GPS  transmitters  in  a  ground- 
based  planar  array.  The  pseudo-GPS  transmitters  are  referred  to  as  pseudolites  in 
this  thesis.  The  second  test  environment  is  required  to  evaluate  the  performance  of 
the  GPS  position  determination  algorithms  developed  in  this  thesis  in  relation  to 
the  performance  achieved  by  conventional  iterative  GPS  algorithms  in  other  than 
typical  GPS  satellite  geometries  where  current  iterative  algorithms  tend  to  perform 
quite  well.  The  test  environment  using  pseudolites  in  ground-based  planar  arrays 
is  indicative  of  the  geometries  that  can  be  encountered  in  a  test  range  such  as  the 
SARS  [25],  [20]  where  areas  of  high  GDOP  or  poor  geometry  are  typical.  These  high 
GDOP  geometrical  conditions  are  not  very  favorable  to  the  conventional  iterative 
GPS  algorithms  and  restrictions  on  flight  profiles  and  ground  array  patterns  must 
be  imposed  to  accommodate  the  deficiencies  in  the  iterative  algorithms  [20].  The 
nonlinear  closed-form  algorithms  are  expected  to  be  more  accurate  in  regions  of  high 
GDOP  than  the  conventional  iterative  algorithms  [2]. 
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The  results  associated  with  the  first  test  environment  are  divided  into  four 
parts,  each  part  related  with  a  different  approach  used  in  obtaining  the  GPS  solu¬ 
tion  estimate.  The  four  parts  include  the  Closed-form  algorithm,  Maximum  likelihood 
algorithm,  Kalman  update  algorithm,  and  Conventional  Iterative  algorithms.  For  the 
second  test  environment,  only  the  results  associated  with  the  Closed-form  algorithm, 
Kalman  update  algorithm,  and  the  Iterative  Least  Squares  (ILS)  algorithm  are  pre¬ 
sented.  All  the  results  presented  are  the  cumulative  representation  of  5000  Monte 
Carlo  runs.  In  order  to  provide  an  unbiased  comparison  basis  between  the  results 
using  the  different  algorithms,  the  Gaussian  pseudorange  noise  for  each  satellite  is 
maintained  the  same  between  the  different  approaches  for  any  given  Monte  Carlo 
run.  The  pseudorange  measurement  noise  for  each  visible  satellite  is  randomly  gener¬ 
ated  as  an  array  with  5000  realizations,  one  for  each  Monte  Carlo  run.  The  simulated 
noise  array  has  a  zero-mean  and  a  standard  deviation  of  cr  =  100  meters  which  is  in 
a  typical  range  but  is  arbitrarily  selected  for  the  purpose  of  the  experimental  trials. 

It  is  not  deemed  necessary  to  carry  out  specific  study  measurements  to  confirm 
the  validity  of  the  selected  standard  deviation  of  cr  =  100  meters  because  the  ex¬ 
perimental  analysis  is  based  strictly  on  relative  performance  of  different  algorithms 
under  identical  noise  conditions.  It  is  nonetheless  essential  that  the  Signal  -  to  -  Noise 
Ratio  (SNR)  used  experimentally  be  in  a  typical  range  since  there  is  a  threshold  level 
below  which  using  stochastic  estimation  theory,  where  an  inverse  covariance  matrix 
is  used  as  a  weighting  matrix  to  obtain  a  minimum  variance  stochastic  estimate,  is 
of  little  or  no  benefit.  Similarly,  if  an  excessively  small  SNR  is  considered  for  exper¬ 
imental  trials,  it  can  lead  to  false  conclusions  because  the  stochastic  approach  will 
outperform  the  deterministic  approach  in  these  unrealistic  conditions.  This  concept 
has  been  demonstrated  by  Brown  for  system  identification  applications  whereby  the 
stochastic  weighted  least  squares  was  found  to  weight  the  estimates  incorrectly  when 
the  SNR  was  above  such  a  threshold  [3].  This  thesis  does  not  attempt  to  determine  a 
SNR  threshold  applicable  specifically  to  GPS  and  of  more  importance  to  determine 
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if  the  current  SNR  encountered  in  typical  GPS  applications  is  above  or  below  this 
threshold.  It  must  be  noted  that  if  the  threshold  is  at  a  SNR  level  much  lower  than 
that  encountered  by  GPS,  it  is  possible  that  there  is  no  benefit  in  applying  stochastic 
estimation  theory. 

Both  GDOP  and  the  condition  number  of  the  the  regressor  matrix  (H),  which 
is  abbreviated  as  cond{H),  are  considered  in  this  thesis  for  the  evaluation  of  how  good 
the  geometrical  conditions  in  the  test  scenarios  are  relative  to  each  other.  GDOP  is 
considered  because  it  is  the  conventional,  most  widely  accepted,  and  most  familiar 
measure  of  geometrical  strength  for  GPS.  The  cond{H)  is  also  being  considered 
because  it  was  shown  that  cond{H)  is  a  better  measure  of  geometrical  strength  since, 
unlike  the  GDOP,  it  has  no  dependence  on  the  satellite  avciilability  [20].  Cond(H)  is 
a  measure  of  strictly  the  geometry;  furthermore,  it  reflects  a  worst  case  upper  limit 
on  error  amplification  that  can  be  produced  by  the  satellite  geometry.  The  condition 
number  of  a  matrix  is  defined  as  the  ratio  of  the  largest  singular  value  to  the  smallest 
singular  value.  In  the  context  of  the  GPS  regressor  matrix  (H)  the  condition  number 
is  calculated  as  follows:  _ 


cond(H)  = 


\\  X^iniHTHY 


where  A(-)  is  used  to  represent  the  eigenvalues  of  a  matrix. 


4.2  Data  Set  Generation,  Simulated  GPS  Test  Environment 

Spilker  [28]  shows  that,  with  a  five  degree  elevation  angle,  at  least  seven  satel¬ 
lites  are  in  view  80  percent  of  the  time  and  there  are  always  at  least  five  satellites  in 
view.  In  terms  of  satellite  availability,  the  worst  case  scenario  occurs  at  latitudes  in 
the  range  of  35  to  55  degrees.  In  this  latitude  range  there  are  at  most  six  satellites 
available  20  percent  of  the  time  using  a  five  degree  elevation  angle.  Even  if  the  worst 
case  conditions  with  a  ten  degree  elevation  angle  were  considered,  in  the  35  to  55 
degree  latitude  range,  only  four  satellites  are  available  less  than  0.5  percent  of  the 
time.  A  ten  degree  elevation  angle  is  extremely  conservative  since  current  technology 
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allows  for  better  than  five  degree  elevation  angles  at  0  altitude.  The  elevation  angle 
can  decrease  to  below  0  degrees  with  increasing  altitude.  Since  satellite  availability 
is  not  dependent  on  user  position  longitude,  selecting  a  single  user  position  in  the  35 
to  55  degree  latitude  range  and  assuming  an  elevation  angle  of  10  degrees  will  allow 
for  simulating  GPS  data  that  is  realistically  indicative  of  worst  case  conditions.  It 
is  important  that  experimental  scenarios  be  diverse  in  addition  to  being  realistic; 
hence,  using  a  fixed  user  position  and  using  the  satellite  availability  at  different  GPS 
times  in  a  24  hour  period  allows  for  a  wide  range  of  diversity. 

An  arbitrary  user  position  in  the  35  to  55  degree  latitude  range  over  the  con¬ 
tinental  United  States,  40°  N  latitude,  105°  W  longitude,  at  an  altitude  of  300  m 
was  selected.  The  geographic  coordinates  are  converted  to  ECEF  coordinates  and 
used  to  generate  the  experimental  data  sets  with  the  Satellite  Navigation  Toolbox 
for  Matlab  [10].  Data  sets  are  generated  at  one  hour  intervals  and  12  scenarios  which 
showed  diversity  in  satellite  availability  and  satellite  geometry  were  selected  to  form 
the  experimental  data  set.  Satellite  availability  ranges  from  five  to  nine,  which  is 
indicative  of  real  life  scenarios.  For  a  satellite  availability  of  six,  seven  and  eight 
where  more  than  one  scenario  was  selected,  the  different  scenarios  were  selected  to 
obtain  the  greatest  diversity  in  satellite  geometries. 

To  give  an  idea  of  the  satellite  geometry  for  each  of  the  twelve  selected  satellite 
scenarios,  a  skyplot  is  generated  for  each  scenario.  The  skyplot  shows  the  relative 
position  of  the  available  satellites  as  seen  at  the  user  position  by  looking  straight  up 
at  the  sky.  The  center  of  the  plot  represents  zenith,  moving  out  from  the  center  of 
the  plot  represents  decreasing  elevation  angles,  and  the  outer  outline  is  the  horizon. 
The  skyplots  for  the  12  experimental  scenarios  are  included  at  Figures  4.6  to  4.17. 

4-3  Simulation  of  the  Ground-Based  Pseudolite  Test  Environment 

The  three  ground-based  pseudolite  test  scenarios  are  produced  based  on  opti¬ 
mized  ground-based  planax  arrays  as  presented  by  McKay  [20].  An  optimum  ground 
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planar  array  of  n  transmitters  has  (n  —  1)  transmitters  uniformly  spread  over  the 
circumference  of  a  circle  of  very  large  radius  ideally  at  the  horizon  and  one  trans¬ 
mitter  at  the  center  of  the  circle.  The  optimum  receiver  location  is  directly  above 
the  transmitter  at  the  center  of  the  circle. 

The  following  three  guidelines  were  used  in  selecting  the  ground  planar  arrays 
for  the  test  environment: 

1.  The  optimum  geometry  for  a  planar  array  based  on  GDOP  as  described  in  [20] 
must  be  maintained  to  avoid  geometries  that  may  potentially  be  very  unfavor¬ 
able  to  the  conventional  iterative  algorithms. 

2.  Maintaining  the  optimum  geometry  concept,  the  three  scenarios  represent  a 
best  possible  condition,  a  degraded  condition  and  a  failure  condition. 

3.  Without  adversely  effecting  GDOP,  sufficient  transmitters  are  used  in  the  pla¬ 
nar  array  to  allow  the  use  of  the  ergodic  assumption  invoked  in  Equation  (3.20) 
to  estimate  the  equation  error  statistics. 

The  test  environment  includes  36  transmitters  uniformly  spread  over  the  cir¬ 
cumference  of  a  circle  of  radius  10000  meters  centered  at  the  origin.  It  is  recognized 
that  a  radius  of  10000  meters  is  restrictive  but  it  is  used  to  reflect  potential  size 
restrictions  due  to  terrain  availability  and/or  line  of  sight  visibility  problems.  The 
user  is  maintained  stationary  directly  above  the  origin  at  an  altitude  of  10000  meters 
(i.e.  user  coordinates  {0,  0,  10000})  for  the  three  scenarios.  To  avoid  confusion  with 
the  scenarios  presented  for  the  first  test  environment,  the  three  scenarios  are  labeled 
Scenarios  A,  B,  and  C.  The  following  differentiates  each  of  the  three  scenarios: 

•  Scenario  A.  The  ground  array  includes  just  the  36  GPS  transmitters  in  the 
outer  circular  pattern. 

•  Scenario  B.  The  ground  array  is  the  same  as  for  Scenario  A  but  an  additional 
transmitter  is  included  at  the  origin  directly  below  the  user  receiver  position. 
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This  scenario  represents  the  lowest  achievable  GDOP  for  the  given  number 
of  receivers  and  ground  area  covered,  or  in  other  words,  the  most  favorable 
conditions  for  the  conventional  ILS  algorithm.  The  planar  array  depicted  in 
Scenario  B  is  shown  in  Figure  4.1. 

•  Scenario  C.  The  ground  array  is  the  same  as  for  Scenario  A,  but  an  additional 
transmitter  is  included  400  meters  away  from  the  origin;  hence,  below  but  not 
directly  below  the  user  receiver  position,  representing  a  degraded  condition  for 
the  conventional  ILS  algorithm. 


Figure  4.1  Scenario  B,  Optimum  Ground-Based  Planar  Array 

4.4  Experimental  Data  Description 

The  following  provides  a  description  of  the  data  collected  using  the  experi¬ 
mental  GPS  position  determination  algorithms.  The  descriptions  are  intended  to 
facilitate  interpretation  of  the  experimental  results.  They  are  presented  in  the  same 
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order  that  they  appear  in  the  tabulated  results  presented  in  Table  4.1  through  Ta¬ 
ble  4.12: 

•  Estimated  a.  This  is  an  estimation  of  the  standard  deviation  of  the  applied 
pseudorange  measurement  noise.  A  data-driven  estimate  of  a  is  calculated  on 
each  Monte  Carlo  run  and  the  result  shown  is  the  average  from  all  the  runs. 
Only  the  closed-form  algorithm  has  the  capability  to  perform  this  calculation. 

•  GDOP.  To  demonstrate  the  diversity  in  geometry,  the  GDOP  for  each  exper¬ 
imental  scenario  is  shown.  GDOP  is  a  measure  of  geometry  associated  with 
the  conventional  iterative  GPS  position  determination  algorithms;  hence,  it  is 
only  presented  in  the  results  for  the  ILS  and  IBF  algorithms.  The  GDOP  is 
calculated  using,  GDOP  =  ^Jtrace{{H'^H)~'^)  where  H  is  the  matrix  of  direc¬ 
tion  cosines  with  ones  populating  the  last  column,  as  defined  in  Section  2.2.2 
of  this  thesis.  The  GDOP  for  the  ILS  algorithm  is  calculated  using  the  tall 
(n  X  4)  AT  matrix  which  reflects  the  n  available  satellites.  The  GDOP  for  the 
IBF  algorithm  is  calculated  using  the  (4  X  4)  if  matrix  which  reflects  the  set 
of  four  satellites  amongst  the  available  satellites,  selected  to  yield  the  lowest 
value  of  GDOP.  In  order  to  obtain  the  lowest  achievable  GDOP,  all  the  possible 
combinations  of  four  satellite  sets  were  tested  for  each  scenario  and  the  best 
combination  was  selected. 

•  Condition  {H).  This  represents  the  condition  number  of  the  regressor  H  used 
in  calculating  the  unknown  GPS  parameters.  For  the  closed-form  algorithm, 
the  effective  regressor  is  considered  in  the  calculation  as,  R~^H^  where 
represents  the  covariance  of  the  measurement  noise  vector  without  the  scalar 
premultiplier  as  defined  in  Equation  (3.11).  Rn  is  used  as  the  noise  strength 
weighting  to  obtain  the  closed-form  estimate  as  defined  in  Equation  (3.17). 
The  calculated  condition  number  is  actually  the  condition  number  of  yjR~^ 
where  the  Cholesky  decomposition  of  R~^  is  used  to  calculate  the  matrix  square 
root.  For  the  iterative  algorithms,  the  regressor  H  is  the  conventional  matrix 
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of  direction  cosines  with  ones  populating  the  last  column.  There  is  no  regres¬ 
sor  associated  with  the  maximum  likelihood  approach  or  the  Kalman  update 
approach. 

•  Error  (x,y,z,bias).  Errors  are  shown  for  user  position  ECEF  coordinates 
(uxtUy^u^)  and  for  the  range-equivalent  user  clock  bias  (6).  For  each  parame¬ 
ter,  the  error  is  the  difference  between  the  sample  average  estimate  produced 
by  the  respective  position  determination  algorithm  and  the  true  value  of  the 
parameter.  The  parameter  sample  average  estimates,  produced  by  the  algo¬ 
rithms,  are  obtained  from  averaging  the  paxameter  estimates  from  the  5000 
Monte  Carlo  runs. 

•  Miss  distance.  The  miss  distance,  or  root  mean  square  (rms)  position  error,  is 
a  measure  of  the  positioning  accuracy  and  is  calculated  from  the  user  position 
ECEF  coordinates  sample  average  estimation  errors.  This  measure  reflects  how 
far  the  estimated  user  position  is  from  the  true  user  position.  The  miss  distance 
is  calculated  using,  md  =  yJerror{xY  -f-  error[yy  -|-  error{zY.  The  estimate 
of  the  range-equivalent  user  clock  bias  has  no  effect  on  the  miss  distance  which 
is  a  function  of  strictly  the  position  parameter  estimation  errors. 

•  Estimated  std  (x,y,z,bias).  The  average  estimated  standard  deviation  is  a 
calculation  associated  with  the  estimate  of  the  GPS  solution  itself.  An  estimate 
of  standard  deviation,  for  each  parameter,  is  calculated  on  each  Monte  Carlo 
run  and  the  result  shown  is  the  average  from  all  the  runs.  This  calculation  is 
related  to  the  estimate  of  a  and  is  performed  within  the  closed-form  algorithm 
only,  since  it  is  the  only  one  with  the  capability.  The  closed-form  algorithm 
is  actually  producing  an  estimate  of  the  covariance  matrix  which  is  a  fully 
populated  matrix  due  to  some  unavoidable  yet  small  and  unknown  level  of  cross 
correlation;  however,  the  off-diagonal  elements  will  be  small.  The  standard 
deviation  for  each  of  the  four  estimated  GPS  parameters  is  obtained  by  taking 
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the  square  root  of  the  diagonal  elements  in  the  predicted  or  estimated  error 
covariance  matrix. 

•  Experimental  std  {x,y,z,bias).  The  experimental  standard  deviation  corre¬ 
sponds  to  the  actual  standard  deviation  of  the  set  of  estimates  collected  over 
all  the  Monte  Carlo  runs.  This  data  is  essentially  the  most  reliable  and  most 
critical  measure  of  relative  performance  of  the  algorithms.  The  experimental 
standard  deviation  is  a  marker  of  how  good  the  experimentally  obtained  esti¬ 
mates  are  in  terms  of  spread,  which  is  a  direct  indication  of  how  sensitive  the 
algorithm  is  to  pseudorange  measurement  noise.  The  standard  deviation  is  a 
more  robust  performance  indicator  than  the  parameter  estimation  errors  or  the 
miss  distance  since,  the  covariance  of  a  random  variable  is  readily  recovered 
using  a  finite  number  of  realizations  and  is  very  consistent.  On  the  other  hand, 
the  mean  is  much  more  difficult  to  recover,  which  raises  some  concern  over  the 
results  for  the  parameter  estimation  errors  and  the  miss  distance,  despite  the 
5000  Monte  Carlo  runs.  The  experimental  standard  deviation  also  serves  as  a 
baseline  to  compare  how  well  the  closed-form  algorithm  predicts  its  parameter 
error  standard  deviations. 

•  Expected  std  (x,  j/,  z,  bias).  The  expected  standard  deviation  is  shown  to  ver¬ 
ify  how  well  the  standard  deviation  of  an  estimate  can  be  calculated  provided 
that  the  standard  deviation  of  the  pseudorange  noise  (cr)  is  known.  The  true 
value  of  cr  is  used  in  the  calculation  even  though  in  reality  the  true  value  of  a 
could  never  be  known.  This  measure  is  an  interesting  one  since  by  comparing 
expected  values  to  the  experimentally  obtained  values,  the  algorithm  used  in 
the  calculation  can  be  validated.  Expected  standard  deviation  is  calculated 
for  all  the  approaches  with  the  exception  of  the  maximum  likelihood  approach 
which  does  not  have  the  capability  to  produce  an  expected  standard  deviation. 
The  only  difference  between  the  expected  standard  deviation  and  the  estimated 
standard  deviation  is  that  the  latter  makes  use  of  a  data-driven  estimate  of 
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the  pseudorange  measurement  noise  variance  which  only  the  closed-form 
algorithm  has  the  capability  of  producing. 

4.5  The  Monte  Carlo  Runs 

The  experimental  trials  consisted  of  5000  Monte  Carlo  runs  for  each  test  sce¬ 
nario.  Although  5000  Monte  Carlo  trials  may  appear  over- abundant,  it  was  observed 
experimentally  that  an  extremely  large  number  of  runs  is  required  in  order  to  recover 
unbiased  means  of  the  estimated  parameters.  Using  experimental  noise  arrays  with 
zero-mean  and  a  standard  deviation  of  <j  =  100  meters,  the  experimental  covariance 
of  the  estimated  parameters  was  found  to  converge  to  a  steady  value  after  a  few 
hundred  runs  but  it  was  not  possible  to  determine  the  number  of  runs  that  would 
be  required  for  the  mean  of  the  estimated  parameters  to  converge  to  a  steady  state 
value.  Under  the  experimental  environments  used  for  this  thesis  it  is  anticipated  that 
over  10000  runs  would  have  been  required.  Running  the  same  experiment  of  5000 
Monte  Carlo  runs  several  times,  using  a  different  zero-mean,  a  =  100  m  simulated 
noise  array  each  time,  produced  virtually  identical  covariance  results  but  some  small, 
yet  significant,  variations  in  the  means  were  still  being  observed. 

4.6  Experimental  Observations,  Simulated  GPS  Test  Environment 

The  four  sections  that  follow  discuss  results  associated  with  the  simulated  GPS 
test  environment.  The  results  are  presented  in  Table  4.1  through  Table  4.12.  The 
results  are  tabulated  by  scenario  and  each  table  is  accompanied  by  a  skyplot  that 
reveals  the  satellite  availability  and  geometry  of  each  scenario. 

4.7  Conventional  Iterative  Algorithms 

The  IBF  and  ILS  algorithms  are  basically  the  conventional  GPS  position  es¬ 
timation  algorithms  currently  being  used.  The  IBF  algorithm  uses  only  a  set  of 
four  satellite  pseudorange  measurements  chosen  from  the  available  satellites  based 
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on  best  GDOP.  This  approach  yields  a  four  by  four  regressor  which  must  be  inverted 
to  obtain  a  solution.  The  ILS  algorithm  is  essentially  the  same  but  all  n  available 
pseudorange  measurements  are  used  to  obtain  the  GPS  solution.  The  regressor  is  a 
(n  X  4)  tall  matrix  so  the  generalized  inverse  is  used  resulting  in  a  least  squares  so¬ 
lution.  The  results  for  these  two  conventional  approaches  are  presented  to  provide  a 
comparison  baseline.  The  ultimate  aim  is  to  improve  on  the  GPS  solutions  produced 
by  the  two  conventional  iterative  algorithms. 

4.7.1  Analysis  of  Results.  The  results  for  both  the  IBF  and  ILS  algorithms 
are  presented  to  provide  a  comparison  baseline  for  the  results  obtained  by  the  closed- 
form,  maximum  likelihood,  and  Kalman  update  algorithms.  Direct  comparison  under 
varying  scenarios  reveals  how  well  the  new  algorithms  are  performing  in  relation  to 
the  conventional  iterative  algorithms. 

Under  some  scenarios,  the  IBF  algorithm  yielded  GPS  solution  estimates  with 
smaller  mean  errors  while  the  ILS  algorithm  yielded  smaller  mean  errors  under  other 
scenarios.  There  was  no  strongly  apparent  geometry  relationship  that  would  indicate 
which  of  the  two  iterative  algorithms  yielded  the  smaller  parameter  estimation  mean 
errors  for  any  given  scenario.  For  both  iterative  algorithms,  a  loose  relationship 
was  observed  between  the  GDOP  and  the  miss  distance,  in  which  the  lower  GDOP 
scenarios  tended  to  yield  smaller  miss  distances.  These  relationships  are  shown  in 
Figure  4.2  which  also  shows  that  the  relationship  extended  to  the  condition  number 
of  H.  It  was  expected  that  a  stronger  relationship  between  GDOP  and  miss  distance 
would  have  existed.  The  lack  of  the  strong  relationship  in  itself  is  an  indicator  of  the 
weakness  in  using  the  parameter  mean  estimation  errors  as  sole  performance  indica¬ 
tors.  It  is  not  possible  to  base  comparisons  conclusively  on  strictly  the  miss  distance 
since  the  results  obtained  show  small  variations  in  the  miss  distances  between  the 
two  iterative  approaches  for  all  the  scenarios.  These  small  differences  in  miss  dis¬ 
tances  may  also  be  partly  due  to  the  biases  in  the  mean  errors  caused  by  insufficient 
Monte  Carlo  runs.  It  is  safer  to  draw  comparisons  on  the  more  robust  experimental 
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standard  deviations  which  tend  to  be  much  more  consistent.  The  relationship  be¬ 
tween  the  standard  deviation  and  GDOP  or  cond{H)  is  extremely  evident,  as  can  be 
seen  in  Figure  4.3  which  is  a  plot  of  the  average  standard  deviation  from  the  three 
position  parameters.  Regardless  of  which  of  the  two  iterative  algorithms  provided 
the  smallest  miss  distance,  the  ILS  algorithm  always  provided  the  smallest  standard 
deviations.  This  implies  that  for  any  given  scenario,  there  is  less  variation  in  the 
results  yielded  by  the  ILS  algorithm,  or  in  other  words,  the  algorithm  is  less  sensitive 
to  the  pseudorange  noise.  For  the  four  estimated  parameters,  in  all  the  scenarios, 
the  standard  deviations  produced  by  the  ILS  algorithm  varied  between  55.56  and 
289.14  meters.  This  will  be  used  as  the  comparison  baseline  for  the  GPS  position 
determination  algorithms  developed  in  this  thesis. 


Iterative  Least  Squares  Algorithm  Iterative  Best  Four  Algorithm 


GDOP  or  cond(H)  GDOP  or  cond(H) 

Figure  4.2  Geometry  Effects  on  Miss  Distance 

The  iterative  algorithms  do  not  have  the  capability  of  predicting  the  estimation 
error  covariance;  however,  when  the  true  value  of  a  is  provided  to  either  of  the 
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conventional  iterative  algorithms,  it  is  possible  to  calculate  the  expected  covariance. 
This  portion  of  the  experimental  trials  is  being  included  for  academic  discussions 
only  since  there  is  no  way  of  knowing  what  the  true  value  of  (T  is;  furthermore,  the 
iterative  algorithms  do  not  have  the  capability  of  producing  a  data-driven  estimate 
of  <7.  Since  the  measurement  noise  on  all  pseudoranges  is  assumed  independent, 
zero-mean,  with  equal  variance  «7^,  the  covariance  of  the  conventional  iterative  GPS 
solution  is  a  function  of  a  given  by,  cov(X)  —  [28].  It  follows  that  the 

expected  standard  deviation  of  the  GPS  solution  parameters  can  be  calculated  using 
the  following  expression,  given  in  Matlab  [17]  notation: 

<^uy,  (Tu,-,  <7fc]  =  \ldiag{cov{X)) 

=  asJdiagiiH^H)-^) 

The  calculated  standard  deviations  of  the  GPS  solution  parameter  estimates  were 
remarkably  close  to  the  values  obtained  experimentally.  For  both  iterative  algo¬ 
rithms,  the  expected  standard  deviations  calculated  where  within  2.25  percent  of  the 
experimentally  obtained  values. 

4-8  Closed-Form  Algorithm 

The  closed-form  algorithm  is  based  on  the  approach  developed  in  Section  3.1 
of  this  thesis.  This  approach  provides  a  closed-form  position  estimate  solution  to  a 
stochastically  modeled  set  of  n  satellite  pseudorange  equations.  The  calculated  user 
position  estimate  is  not  dependent  on  the  equation  error  covariance,  Py,  which  can 
also  be  interpreted  as  the  noise  strength;  however,  the  noise  strength  is  required  in 
the  derivation  of  the  estimate.  This  approach  also  provides  an  estimate  of  the  co- 
variance  of  the  estimation  error  using  the  return  dilferences  or  residuals.  The  return 
difference  vector,  as  defined  in  Section  3.1.4,  is  used  to  obtain  a  data-driven  estimate 
of  the  pseudorange  measurement  noise  variance  (cr^),  from  which  the  estimation  error 
covariance  is  readily  calculated. 
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Iterative  Least  Squares  Algorithm 


Iterative  Best  Four  Algorithm 
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Figure  4.3  Geometry  Effects  on  Standard  Deviation 


4.8.1  Limitations.  Some  known  limitations  associated  with  the  closed- 
form  algorithm  must  be  emphasized  to  avoid  misinterpretation  of  the  experimental 
results.  The  following  are  the  limitations  that  must  be  kept  in  mind  in  analyzing 
the  experimental  results: 


1.  n  satellite  pseudorange  measurements  are  required  to  form  a  linear  regression 
of  n  —  1  equations.  Since  there  are  four  unknown  parameters,  the  user  position 
ECEF  coordinates  and  the  range-equivalent  user  clock  bias,  a  minimum  of  five 
pseudorange  measurements  are  required. 

2.  When  exactly  five  pseudorange  measurements  are  available  (i.e.  n  =  5),  the 
linear  regression  takes  the  form  of  four  equations  in  four  unknowns.  In  this  sit¬ 
uation,  return  differences  of  essentially  zero  can  be  expected,  thereby  rendering 
the  covariance  estimation  capability  of  the  algorithm  invalid. 
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3.  In  addition  to  the  previous  items,  n  must  be  large  to  determine  a  experimen¬ 
tally.  It  is  not  known  how  large  n  must  be  before  a  reliable  estimate  of  cr  can  be 
extracted  from  the  return  difference,  but  an  unconfirmed  rule  of  thumb  from 
statistics  suggests  that  =  2iV^,  where  N  is  the  number  of  variables.  For 
GPS,  =  4  which  implies  an  n^m  of  at  least  32.  The  value  of  serves 
only  as  a  guideline  because  it  is  desirable  to  have  n  as  large  as  possible  since 
the  ergodic  assumption  invoked  in  Equation  (3.20)  becomes  more  valid  as  the 
number  of  available  satellites  increases.  The  impact  of  this  limitation  can  only 
be  assessed  through  thorough  experimentation. 

4-8.2  Analysis  of  Results.  The  experimental  results  indicate  that  the 
closed-form  algorithm  developed  in  this  thesis  is  extremely  sensitive  to  noise.  The 
sensitivity  to  pseudorange  noise  is  reflected  in  the  extremely  large  experimental 
standard  deviations  which  ranged  from  1078.08  m  to  23864.01  m  for  the  position 
parameter  estimates  as  compared  to  ILS  algorithm’s  baseline  results  of  55.56  m  to 
289.14  m.  The  standard  deviations  for  the  range-equivalent  user  clock  bias  esti¬ 
mates  were  significantly  larger  ranging  between  19014.42  m  and  111746.08  m.  For 
the  twelve  scenarios  under  observation,  the  miss  distance  ranged  between  4.93  m 
and  729.17  m  compared  to  the  baseline  ILS  results  of  0.76  m  to  7.84  m.  The  re¬ 
sults  loosely  reflect  tendencies  that  better  results  (i.e.  lower  mean  errors  and  miss 
distance)  are  obtained  when  the  condition  number  of  the  regressor  H  is  lower.  The 
relationship  between  the  condition  number  of  the  regressor  H  and  miss  distance  is 
plotted  in  Figure  4.4.  In  the  twelve  test  scenarios  the  results  both  in  terms  of  mean 
estimation  errors  and  error  standard  deviations  were  significantly  worse  than  those 
obtained  by  both  conventional  iterative  algorithms. 

The  results  show  that  the  condition  numbers  for  H  yielded  by  the  closed-form 
algorithm  are  all  extremely  large,  ranging  between  220.92  and  2172.53;  however, 
an  in  depth  appreciation  of  the  algorithm  has  not  yet  been  acquired  to  allow  in 
distinguishing  what  could  be  considered  a  good  range  for  the  condition  numbers. 
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Figure  4.4  Closed- Form  Algorithm  Performance  versus  Condition  Number 

Consequently,  it  is  unlikely  that  this  algorithm  can  provide  a  GPS  solution  with 
small  errors.  The  ill  conditioning  of  the  regressor  is  largely  due  to  the  last  column 
which  is  made  up  of  the  difference  between  the  pseudoranges.  The  first  three  columns 
are  made  up  of  differences  in  the  three  satellite  position  coordinates  (uj,,  Uy,  u^), 
respectively,  which  tend  to  produce  much  larger  differences  than  the  pseudorange 
differences.  Based  on  the  formulation  of  the  H  matrix,  the  lower  condition  numbers 
are  achieved  when  the  spread  in  the  pseudorange  measurements  of  the  satellites 
in  view  is  large.  If  all  pseudoranges  are  equal,  which  is  not  a  physically  possible 
scenario  in  the  NAVSTAR  GPS  satellite  constellation,  the  last  column  of  the  H 
matrix  would  be  all  zeros,  yielding  a  condition  number  of  infinity.  The  poor  scaling 
due  to  the  last  column  of  the  H  matrix  manifests  itself  as  extremely  large  errors 
in  the  range-equivalent  user  clock  bias  where  the  observed  errors  ranged  between 
2.04  m  and  2500.71  m  compared  to  baseline  results  of  0.03  m  to  5.52  m.  The  errors 
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in  the  estimation  of  the  range-equivalent  user  clock  bias  were  generally  in  the  range 
of  one  order  of  magnitude  larger  than  the  errors  observed  for  the  three  coordinates  of 
the  user  position  estimate.  The  large  errors  in  the  clock  bias  estimates  do  not  affect 
the  miss  distance  which  is  strictly  a  function  of  the  error  in  the  estimated  position. 

An  additional  feature  of  the  closed-form  algorithm  is  its  ability  to  provide 
a  data-driven  estimate  of  the  covariance  of  the  GPS  solution  estimate.  This  is 
achieved  by  using  the  return  difference,  or  residuals,  to  estimate  the  variance  of  the 
pseudorange  measurement  noise  (cr^).  It  is  recognized  that  when  only  five  satellite 
pseudorange  measurements  are  available,  it  is  not  possible  to  estimate  the  error 
covariance  since  the  the  return  difference  will  be  virtually  zero.  The  results  show 
that,  with  more  than  five  satellites,  an  estimate  of  <7,  the  standard  deviation  of  the 
experimental  pseudorange  measurement  noise,  is  obtained.  The  estimates  of  cr  are 
not  reliable  due  to  the  small  number  of  satellites  available.  As  satellite  availability 
increases,  the  estimate  of  a  improves  accordingly.  The  experimental  results  show 
that,  with  six  satellites  in  view,  the  estimated  a  ranged  between  23.51  m  and  39.08  m. 
With  seven  satellites,  the  estimated  a  ranged  between  52.87  m  and  56.26  m,  showing 
some  improvement.  With  eight  satellites,  the  estimated  a  ranged  between  62.95  m 
and  64.01  m,  showing  some  further  improvement.  Comparison  of  the  estimated 
standard  deviations  obtained  to  the  actual  cr  of  100  m  used  in  the  generation  of 
the  pseudorange  measurement  noise  clearly  indicates  that  a  much  larger  number  of 
satellites  must  be  available  before  a  reliable  estimate  of  cr  can  be  extracted  from 
the  return  difference.  This  can  present  a  serious  obstacle  since  there  are  very  rarely 
more  than  nine  satellites  in  view  from  the  NAVSTAR  GPS  constellation.  This 
situation  can  change  if  consideration  is  given  to  GLONASS  satellite  augmentation 
or  pseudolite  augmentation. 

The  estimates  of  the  error  covariance  yielded  by  the  closed-form  algorithm  were 
poor  as  a  direct  consequence  of  the  algorithm’s  inability  to  estimate  a  accurately. 
There  are  large  differences  between  the  actual  experimental  standard  deviations  and 
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the  estimated  standard  deviations;  however,  this  was  not  due  to  any  deficiency  in 
the  algorithm  used  in  predicting  the  estimation  error  covariance.  Although  there 
would  not  be  any  way  of  knowing  the  true  a  in  the  real  world,  in  order  to  validate 
the  algorithm  used  in  estimating  the  error  covariance,  the  actual  <t  of  100  m  was 
provided  for  use  in  the  calculations  of  the  expected  standard  deviation.  The  expected 
standard  deviation  calculations  yielded  numbers  extremely  close  to  those  observed 
experimentally.  The  expected  standard  deviations  where  all  within  6.5  percent  of 
the  experimentally  obtained  values.  This  is  sufficient  to  validate  that  the  algorithm 
can  be  used  to  estimate  the  standard  deviation  of  its  GPS  solution  estimate,  but  it 
does  not  change  the  fact  that  the  algorithm  is  plagued  with  extremely  large  GPS 
parameter  estimate  mean  errors  with  extremely  large  standard  deviations.  It  must 
be  noted  that  when  the  true  value  of  a  is  provided  to  either  of  the  conventional 
iterative  algorithms,  the  expected  standard  deviations  of  the  GPS  unknown  parameter 
estimates  obtained  were  also  remarkably  close  to  the  values  obtained  experimentally. 

In  order  to  appreciate  the  merits  of  the  closed-form  algorithm  truly,  particularly 
the  ability  to  predict  the  estimation  error  covariance,  it  is  necessary  to  increase  the 
satellite  availability.  Nothing  can  be  done  to  increase  the  actual  satellite  availability 
which  is  directly  dependent  on  the  constellation,  but  by  applying  some  creative  ex¬ 
perimentation,  it  may  be  possible  to  increase  the  effective  satellite  availability.  If  it 
is  assumed  that  the  noise  is  predominantly  receiver  related  and  it  is  slowly  varying, 
the  estimate  a  can  be  generated  every  k  sample  periods  using  a  moving  window  type 
approach,  in  which  only  the  most  recent  k  data  sets  are  used  to  establish  the  noise 
strength.  This  concept  is  actually  more  complex  than  it  appears  because,  although  it 
is  possible  to  assume  that  the  noise  strength  remains  fairly  constant  over  the  window, 
the  vehicle  is  moving  continuously.  This  concept  has  not  yet  been  evaluated  experi¬ 
mentally;  however,  it  proved  to  be  extremely  effective  and  reliable  for  estimating  a 
while  the  user  position  is  maintained  stationary.  The  experimental  trials  involved 
producing  a  concatenation  of  satellite  scenarios  generated  at  one  hour  time  inter- 
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vals.  The  trials  where  run,  incrementing  the  number  of  concatenated  scenarios  to 
determine  the  effective  satellite  availability  required  to  produce  reliable  estimates  of 
a.  Experimental  results  showed  the  estimated  a  improving  with  increasing  effective 
satellite  availability  until  a  steady  state  maximum  was  reached  using  six  scenarios, 
giving  an  effective  satellite  availability  of  42.  These  results  are  shown  graphically  in 
Figure  4.5.  An  estimated  mean  a  of  92.5  m  was  retrieved  over  1000  Monte  Carlo 
runs.  An  experimental  standard  deviation  of  30.0  m  was  observed  on  the  estimates 
of  <T,  which  implies  that  68.3  percent  of  the  time  the  algorithm  produces  an  estimate 
of  <7  between  62.5  m  and  122.5  m.  Although,  the  user  position  is  stationary  here, 
which  is  seldom  the  case  in  reality,  it  can  be  inferred  that  the  concept  of  scenario 
concatenation  can  be  modified  and  applied  to  increase  the  effective  satellite  avail¬ 
ability  even  while  the  user  position  is  not  stationary.  Being  able  to  estimate  cr,  even 
if  only  while  the  receiver  is  stationary,  can  prove  to  be  extremely  useful  for  receiver 
calibration  purposes. 


Effective  Satellite  Availability 

Figure  4.5  Estimation  of  a  Using  Scenario  Concatenation 
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4-8.3  Experimental  Results,  Alternate  Closed-Form  Algorithm.  The  ex¬ 
perimental  results  obtained  with  the  alternate  closed-form  algorithm  developed  in 
Section  3.4  are  discussed.  The  alternate  closed-form  algorithm,  as  its  name  implies, 
is  simply  a  variation  on  the  closed-form  algorithm  which  was  motivated  by  the  de¬ 
sire  to  obtain  a  more  mathematically  balanced  algorithm  that  can  produce  a  solution 
with  only  four  pseudorange  measurements  as  is  the  case  for  deterministic  algorithms. 
The  mathematical  development  showed  that,  due  to  the  singularity  in  the  noise  co- 
variance  matrix,  a  minimum  of  five  pseudorange  measurements  is  required  just  as 
was  the  case  in  the  closed-form  solution.  Given  the  apparent  similarities  between 
both  closed-form  algorithms,  their  relative  performances  were  not  expected  to  be 
much  different.  In  fact,  in  all  12  experimental  scenarios,  the  experimental  results 
yielded  by  both  versions  of  the  alternate  closed-form  algorithm  were  identical  to 
those  obtained  with  the  closed-form  algorithm.  The  two  versions  of  the  alternate 
closed-form  algorithm  refers  to: 

1.  The  solution  in  Equation  (3.71)  based  on  the  (n  x  n)  linear  regression  with  the 
singular  noise  covariance  matrix. 

2.  The  solution  in  Equation  (3.77)  based  on  the  transformed  (n  —  1  x  n  —  1)  linear 
regression  with  the  diagonal  noise  covariance  matrix. 

The  experimental  results  produced  by  the  alternate  closed-form  algorithm  are  not 
presented  explicitly  since  they  are  identical  to  the  closed-form  results;  hence,  the 
closed-form  results  presented  in  Tables  4.1  through  4.12  can  be  interpreted  as  being 
the  results  from  the  alternate  closed-form  algorithm  as  well.  Contrary  to  the  origi¬ 
nal  expectations,  the  alternate  closed-form  algorithm  does  not  offer  any  additional 
benefits  over  the  the  closed-form  algorithm;  hence,  given  its  complex  structure  and 
the  inherited  singularities  in  the  equation  error  covariance  and  regressor  matrices, 
its  use  is  not  recommended. 
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4-9  Maximum  Likelihood 

The  maximum  likelihood  algorithm  was  developed  as  a  supplement  to  the 
closed-form  algorithm  and  is  not  intended  to  provide  a  GPS  solution  on  its  own. 
The  maximum  likelihood  algorithm  is  initiated  when  the  closed-form  algorithm  has 
produced  a  solution  and  the  closed-form  solution  is  used  as  a  starting  point  or  ini¬ 
tialization  for  the  maximum  likelihood  algorithm.  The  estimate  of  the  unknown 
GPS  parameter  vector  and  the  associated  error  covariance  obtained  from  the  closed- 
form  algorithm  is  sufficient  to  describe  the  Gaussian  probability  distribution  func¬ 
tion  entirely.  A  constrained  optimization  problem  is  formulated  by  using  the  last 
pseudorange  equation  at  Equation  (3.3)  which  was  not  fully  employed  in  the  linear 
regression.  The  goal  of  the  optimization  problem  is  to  obtain  a  new  estimate  of  the 
GPS  parameters  that  maximizes  the  probability  subject  to  the  constraint  imposed 
by  the  pseudorange  measurement  equation;  hence,  the  name  maximum  likelihood. 

4.9.1  Limitations.  The  following  are  some  known  limitations  associated 
with  the  maximum  likelihood  algorithm  which  must  be  kept  in  mind  in  analyzing  the 
experimental  results: 

1.  Optimization  algorithms  are  iterative  in  nature,  hence  the  closed-form  appeal 
is  lost. 

2.  The  maximum  likelihood  algorithm  can  provide  an  improved  estimate  of  the 
GPS  unknown  parameter  vector,  but  there  is  no  simple  way  of  predicting  the 
estimation  error  covariance. 

3.  The  pseudorange  measurement  associated  with  the  last  satellite  is  noise  cor¬ 
rupted  but  the  constraint  equation  used  in  the  development  of  the  algorithm 
must  be  a  deterministic  equation.  Disregarding  the  noise  in  the  equation 
will  bias  the  solution  and  can  have  an  adverse  impact  on  the  error  since  the 
solution  becomes  extremely  dependent  on  the  noise  of  that  one  single  pseudo¬ 
range  measurement.  Some  of  this  concern  is  alleviated  by  recalling  that  the  last 
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pseudorange  equation  is  actually  the  one  closest  to  the  average  pseudorange  as 
defined  in  Equation  (3.4).  As  such,  it  is  likely  that  the  last  satellite  is  some¬ 
what  in  the  vicinity  of  zenith  where  lower  levels  of  pseudorange  measurement 
noise  can  be  expected. 

4.9.2  Analysis  of  Results.  The  results  produced  by  the  maximum  likeli¬ 
hood  algorithm  demonstrate  the  necessity  of  using  the  last  GPS  pseudorange  mea¬ 
surement.  In  all  twelve  scenarios,  regardless  of  what  the  GPS  solution  errors  where 
when  the  algorithm  was  initiated,  the  errors  are  reduced  drastically  to  essentially 
the  same  level  as  obtained  with  ILS  method.  The  most  remarkable  difference  is  in 
the  errors  of  the  range-equivalent  user  clock  bias  which  were  extremely  large  in  the 
initial  conditions  provided.  There  is  no  direct  relationship  between  the  closed-form 
results  and  the  maximum  likelihood  results  in  that  the  smaller  initial  errors  do  not 
necessarily  yield  better  results.  The  miss  distances  yielded  ranged  from  0.89  m  to 
8.22  m  and  followed  the  ILS  results  quite  closely  for  the  four  GPS  solution  param¬ 
eters  in  all  twelve  scenarios.  The  miss  distances  yielded  by  the  maximum  likelihood 
algorithm,  when  compared  to  those  yielded  by  the  ILS  algorithm,  are  a  fraction  of 
a  meter  either  above  or  below  and  the  largest  difference  was  0.56  m  observed  in 
Scenario  8. 

The  experimental  standard  deviations  observed  for  the  three  position  estima¬ 
tion  parameters  are  essentially  the  same  as  those  observed  for  the  ILS.  The  per¬ 
formance  in  terms  of  to  the  estimation  of  the  range-equivalent  user  clock  bias  is 
slightly  worse  than  the  baseline  established  by  the  ILS  algorithm  both  in  terms  of 
the  mean  error  and  the  error  covariance.  If  only  the  user  position  estimation  param¬ 
eters,  {ux,Uy,Uz)j  are  considered,  the  maximum  likelihood  algorithm  demonstrates 
performance  equivalent  to  the  baseline.  The  fundamental  limitation  with  the  max¬ 
imum  likelihood  approach  is  its  inability  to  predict  its  estimation  error  covariance; 
furthermore,  unlike  the  iterative  approaches,  even  if  the  standard  deviation  of  the 
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pseudorange  measurement  noise  is  somehow  known,  the  expected  standard  deviation 
of  the  estimation  parameters  cannot  be  calculated. 

The  results  yielded  by  the  maximum  likelihood  algorithm  provide  some  hope 
that  a  stochastic  estimation  approach  to  the  GPS  pseudorange  equations  is  feasi¬ 
ble  and  that  performance  comparable  to  that  currently  being  achieved  by  the  ILS 
algorithm  is  attainable.  However,  the  limitations  associated  with  the  maximum  like¬ 
lihood  algorithm  outweigh  the  advantages  that  motivated  this  research  in  the  first 
place.  Consequently,  another  supplementary  algorithm  that  is  noniterative  and  more 
elegant  in  form  and  that  can  refine  the  solution  produced  by  the  closed-form  algo¬ 
rithm  while  maintaining  the  capability  to  provide  the  estimation  error  covariance  is 
needed,  which  leads  to  the  Kalman  update  algorithm  presented  in  Section  3.3. 

4.10  Kalman  Update 

The  Kalman  update  algorithm  was  developed  as  an  alternate  supplement  to  the 
closed-form  algorithm  and,  as  was  the  case  for  the  maximum  likelihood  algorithm,  it 
is  not  intended  to  provide  a  GPS  solution  on  its  own.  The  Kalman  update  algorithm 
requires  the  closed-form  solution  for  initialization.  The  last  pseudorange  equation  at 
Equation  (3.3)  which  was  not  fully  employed  in  the  linear  regression  is  manipulated 
into  the  form  of  a  new  linear  measurement  and  is  used  to  update  the  solution  yielded 
by  the  closed-form  algorithm  through  an  update  cycle  closely  resembling  that  of  a 
conventional  Kalman  filter. 

4.10.1  Analysis  of  Results.  The  Kalman  update  algorithm  produced  results 
comparable  to  the  baseline  results.  Again,  the  errors  on  the  estimation  parameters 
and  the  miss  distance  obtained  are  extremely  close  to  the  baseline  results.  The  miss 
distances  yielded  by  the  Kalman  Update  algorithm  when  compared  to  those  yielded 
by  the  ILS  algorithm  are  a  fraction  of  a  meter  either  above  or  below  and  the  largest 
difference  was  0.77  m  observed  in  Scenario  8.  Given  the  lack  of  confidence  that  can 
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be  placed  on  just  the  mean  error  results,  the  results  are  considered  equivalent  to  the 
baseline  in  terms  of  parameter  estimation  errors.  For  the  three  position  estimation 
parameters,  the  experimental  standard  deviations  are  also  equivalent  to  the  baseline; 
however,  the  standard  deviation  associated  with  the  error  on  the  range-equivalent 
user  clock  bias  is  noticeably  larger  than  the  baseline.  The  standard  deviation  on  the 
user  clock  bias  was  consistently  marginally  larger  than  those  obtained  with  the  max¬ 
imum  likelihood  algorithm.  Unlike  the  maximum  likelihood  algorithm,  the  Kalman 
update  algorithm  has  the  capability  of  predicting  its  estimation  error  covariance,  and 
its  performance  is  reasonably  good.  The  calculated  standard  deviations  for  the  three 
position  estimation  parameters  where  all  within  14.02  percent  of  the  experimentally 
obtained  values.  The  performance  in  producing  the  expected  standard  deviation  on 
the  user  clock  bias  estimation  error  was  poorer  and  differed  from  the  experimentally 
obtained  values  by  as  much  as  44.2  percent. 

Overall,  the  performance  of  the  Kalman  update  algorithm  is  good  and  truly 
deserving  of  further  consideration.  The  only  shortcoming  of  this  algorithm  is  that  it 
produces  user  clock  bias  estimation  error  standard  deviations  larger  than  the  baseline 
and  the  algorithm  is  incapable  of  predicting  the  standard  deviation,  associated  with 
the  user  clock  bias  estimation  error,  accurately.  For  the  most  part,  the  performance 
of  the  Kalman  update  algorithm  is  comparable  to  that  of  the  baseline;  furthermore, 
the  algorithm  retained  all  the  attractive  features  that  motivated  the  development  of 
the  closed-form  algorithm.  Looking  at  the  closed-form  algorithm  as  supplemented 
by  the  Kalman  update  algorithm  as  a  single  GPS  position  determination  algorithm 
has  produced  an  algorithm  with  the  following  attributes: 

1.  The  algorithm  is  closed-form,  hence  it  can  be  used  under  any  geometrical 
conditions  without  the  need  for  externally  provided  initialization. 

2.  The  algorithm  has  the  potential  to  benefit  from  computational  efficiencies  due 
to  its  non-recursive  nature. 
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3.  The  algorithm  has  the  capability  to  produce  an  estimate  of  the  estimation  error 
covariance. 

4.  The  algorithm  can  produce  a  GPS  solution  estimate  without  the  knowledge  of 
the  pseudorange  measurement  noise  standard  deviation  (cr). 

5.  The  performance  under  typical  navigation  scenarios,  using  only  the  NAVSTAR 
GPS  satellite  constellation,  is  for  the  most  part  equivalent  to  the  performance 
achieved  by  the  conventional  ILS  algorithm  used  as  a  baseline. 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  20700  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  9 
Satellites  for  best  GDOP:  5 17  20  23 
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Figure  4.6  Skyplot,  Scenario  1 


Table  4.1  Monte  Carlo  Simulation  Results,  Scenario  1 


Iterative 

Closed 

Maximum 

Kalman 

Least 

Best 

Form 

Likelihood 

Update 

Squares 

Four 

Estimated  a 

61.69 

GDOP 

1.81 

2.45 

Condition  H 

482.80 

5.43 

Msm 

Error  x 

25.50 

0.36 

0.29 

0.21 

Error  y 

94.51 

-0.17 

0.07 

-0.64 

-0.82 

Error  z 

-83.29 

-1.13 

-1.33 

-0.77 

0.38 

Error  bias 

453.44 

1.12 

2.25 

-0.03 

0.02 

Miss  distance 

1.40 

Estimated  std  x 

Estimated  std  y 

2707.30 

Estimated  std  z 

2349.87 

Estimated  std  bias 

12930.87 

Experimental  std  x 

1137.75 

56.14 

56.35 

55.56 

107.56 

Experimental  std  y 

4289.53 

99.39 

101.73 

97.07 

140.85 

Experimental  std  z 

3723.09 

114.77 

140.61 

Experimental  std  bias 

20433.70 

84.56 

95.91 

Expected  std  x 

1166.00 

57.39 

Expected  std  y 

91.82 

Expected  std  z 

105.55 

■KiaWM 

Expected  std  bias 

20959.97 

120.56 

83.93 

95.10 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  38700  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  8 
Satellites  for  best  GDOP:  1  19  20  22 
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Figure  4.7  Skyplot,  Scenario  2 


Table  4.2  Monte  Carlo  Simulation  Results,  Scenario  2 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

64.01 

GDOP 

1.82 

2.09 

Condition  H 

4.95 

3.63 

Error  x 

-0.60 

-0.68 

-0.77 

Error  y 

37.60 

-1.38 

-1.31 

Error  z 

-32.21 

1.27 

1.21 

Error  bias 

185.69 

-0.03 

0.35 

0.30 

-0.07 

Miss  distance 

50.47 

1.98 

1.88 

2.17 

Estimated  std  x 

707.00 

Estimated  std  y 

2637.00 

Estimated  std  z 

2264.56 

Estimated  std  bias 

12554.94 

Experimental  std  x 

1078.08 

64.09 

64.35 

msKlm 

Experimental  std  y 

4011.12 

122.82 

124.25 

Experimental  std  z 

3442.63 

96.59 

Experimental  std  bias 

19014.42 

113.86 

75.56 

Expected  std  x 

1104.56 

60.16 

64.12 

85.13 

Expected  std  y 

4119.82 

126.32 

10^ 

Expected  std  z 

3537.94 

89.79 

■iVJifeM 

Expected  std  bias 

19614.72 

106.06 

75.51 

80.72 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  13500  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  8 
Satellites  for  best  GDOP:  3 10  15  24 
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Figure  4.8  Skyplot,  Scenario  3 


Table  4.3  Monte  Carlo  Simulation  Results,  Scenario  3 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

62.95 

GDOP 

1.85 

2.39 

Condition  H 

1071.20 

4.49 

Error  x 

-23.65 

0.87 

0.86 

0.98 

Error  y 

-92.07 

1.73 

1.69 

1.69 

3.88 

Error  z 

79.90 

-0.35 

-0.32 

-0.36 

-1.58 

Error  bias 

-440.27 

-0.73 

-0.92 

-1.28 

-2.15 

Miss  distance 

124.18 

1.97 

1.92 

1.94 

4.31 

Estimated  std  x 

1291.51 

Estimated  std  y 

4938.60 

Estimated  std  z 

4225.10 

Estimated  std  bias 

23139.21 

Experimental  std  x 

1929.29 

68.56 

68.13 

82.90 

Experimental  std  y 

7374.58 

125.97 

123.81 

161.19 

Experimental  std  z 

6308.96 

88.19 

86.69 

Experimental  std  bias 

34510.60 

109.77 

84.05 

■niwiM 

Expected  std  x 

2051.48 

63.71 

Expected  std  y 

7844.66 

140.00 

122.84 

Expected  std  z 

6711.31 

89.08 

87.20 

120.19 

Expected  std  bias 

36755.24 

86.79 

83.42 

100.47 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  6300  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  7 
Satellites  for  best  GDOP:  10 1 1 15  21 
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Figure  4.9  Skyplot,  Scenario  4 


Table  4.4  Monte  Carlo  Simulation  Results,  Scenario  4 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

54.13 

GDOP 

1.89 

2.34 

Condition  H 

483.67 

4.81 

Error  x 

-30.88 

-1.44 

-1.48 

Error  y 

bees 

0.32 

BSH 

Error  z 

BSdH 

-0.13 

mm 

Error  bias 

-528.90 

-1.00 

-1.64 

-0.06 

Miss  distance 

BIsSH 

1.48 

Estimated  std  x 

Estimated  std  y 

3263.67 

Estimated  std  z 

2827.15 

Estimated  std  bias 

15479.98 

Experimental  std  x 

1619.75 

68.97 

69.09 

Experimental  std  y 

6121.03 

116.08 

117.33 

114.85 

130.93 

Experimental  std  z 

5304.14 

106.57 

104.87 

130.06 

Experimental  std  bias 

29006.92 

105.45 

81.97 

89.33 

Expected  std  x 

1596.27 

65.01 

69.37 

110.06 

Expected  std  y 

6029.36 

130.54 

Expected  std  z 

5222.94 

113.60 

104.46 

132.20 

Expected  std  bias 

28598.00 

96.66 

82.28 

90.70 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  56700  sec 
User  Position:  40  N.  105  W,  300  m 
Satellite  Availability:  N  =  7 
Satellites  for  best  GDOP:  3  4  9  13 
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Figure  4.10  Skyplot,  Scenario  5 
Table  4.5  Monte  Carlo  Simulation  Results,  Scenario  5 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

56.26 

GDOP 

Condition  H 

850.41 

5.70 

5.63 

Error  x 

46.31 

-0.25 

-0.34 

-0.67 

Error  y 

172.78 

2.11 

2.17 

1.93 

3.12 

Error  z 

-146.70 

-0.26 

-0.31 

-0.08 

-1.57 

Error  bias 

■rffekfiM 

-0.43 

-0.13 

-0.75 

-1.90 

Miss  distance 

2.14 

2.20 

1.96 

3.55 

Estimated  std  x 

Estimated  std  y 

5493.08 

Estimated  std  z 

4713.29 

Estimated  std  bias 

25600.41 

Experimental  std  x 

2694.47 

76.40 

76.47 

Experimental  std  y 

9857.39 

134.06 

134.25 

Experimental  std  z 

8461.39 

111.63 

162.96 

Experimental  std  bias 

45921.95 

115.69 

99.73 

126.19 

Expected  std  x 

2665.27 

83.14 

Expected  std  y 

9763.41 

150,57 

134.95 

170.91 

Expected  std  z 

8377.40 

110.55 

108.26 

159.81 

Expected  std  bias 

45502.20 

98.12 

98.79 

125.21 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time;  t  =  85500  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  7 
Satellites  for  best  GDOP:  8121518 
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Figure  4.11  Skyplot,  Scenario  6 


Table  4.6  Monte  Carlo  Simulation  Results,  Scenario  6 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  cr 

55.60 

GDOP 

Kiyn 

kiei 

Condition  H  i 

792.75 

7.01 

Error  x 

-2.58 

-2.32 

-2.32 

-2.33 

nog 

Error  y 

-3.47 

-2.53 

-2.54 

-2.49 

Error  z 

2.36 

1.54 

1.75 

Error  bias 

-2.04 

2.39 

Miss  distance 

4.93 

3.76 

3.77 

3.73 

3.91 

Estimated  std  x 

1148.74 

Estimated  std  y 

4132.45 

Estimated  std  z 

3614.97 

Estimated  std  bias 

19524.84 

Experimental  std  x 

2053.96 

87.66 

87.67 

87.49 

117.39 

Experimental  std  y 

7398.57 

124.98 

Experimental  std  z 

6462.19 

154.15 

154.18 

153.81 

190.61 

Experimental  std  bias 

34910.46 

134.45 

135.37 

Expected  std  x 

2066.11 

91.89 

86.96 

115.34 

Expected  std  y 

7432.62 

127.00 

124.90 

163.45 

Expected  std  z 

6501.87 

■IriililM 

153.29 

191.49 

Expected  std  bias 

35117.32 

122.17 

147.89 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  53100  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  7 
Satellites  for  best  GDOP:  4  7  14  21 
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Figure  4.12  Skyplot,  Scenario  7 


Table  4.7  Monte  Carlo  Siitiulation  Results,  Scenario  7 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

52.87 

GDOP 

2.96 

3.90 

Condition  H 

2179.74 

8.67 

8.66 

Error  x 

0.04 

0.06 

Error  y 

0.22 

0.26 

Error  z 

-48.62 

-1.33 

H&B 

Error  bias 

251.92 

-0.89 

-1.12 

Miss  distance 

74.32 

1.38 

1.34 

1.42 

0.35 

Estimated  std  x 

2621.91 

Estimated  std  y 

9811.00 

Estimated  std  z 

8569.22 

Estimated  std  bias 

45828.57 

Experimental  std  x 

5096.16 

69.62 

69.61 

69.56 

100.28 

Experimental  std  y 

19066.07 

175.83 

198.45 

Experimental  std  z 

16649.71 

168.38 

167.98 

251.23 

Experimental  std  bias 

89017.21 

162.90 

167.95 

194.91 

Expected  std  x 

4958.95 

65.48 

69.45 

101.08 

Expected  std  y 

18556.05 

173.72 

174.60 

199.25 

Expected  std  z 

16207.40 

Expected  std  bias 

86677.93 

146.43 

153.94 

196.20 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  45900  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  6 
Satellites  for  best  GDOP:  4  6 14  22 


Figure  4.13  Skyplot,  Scenario  8 


Table  4.8  Monte  Carlo  Simulation  Results,  Scenario  8 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  cr 

23.51 

bhi 

GDOP 

2.84 

Condition  H 

337.28 

7.45 

itftaiii 

Error  x 

-19.42 

-0.23 

-0.14 

Error  y 

-73.93 

-0.85 

-0.18 

Error  z 

64.35 

1.12 

1.25 

0.73 

Error  bias 

-347.06 

-0.98 

-1.73 

-0.29 

Miss  distance 

99.92 

1.32 

1.53 

0.76 

1.53 

Estimated  std  x 

287.39 

Estimated  std  y 

1094.28 

Estimated  std  z 

944.04 

Estimated  std  bias 

5165.26 

1 

1 

Experimental  std  x 

1246.90 

66.20 

66.71 

65.55 

Experimental  std  y 

4738.02 

216.20 

217.87 

mmoi 

236.33 

Experimental  std  z 

4084.66 

108.34 

110.45 

125.90 

Experimental  std  bias 

22298.07 

144.68 

171.07 

150.91 

Expected  std  x 

1222.25 

67.52 

85.50 

Expected  std  y 

4653.81 

221.57 

236.35 

Expected  std  z 

4014.87 

100.53 

Expected  std  bias 

21967.14 

126.23 

136.69 

151.20 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  71 100  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  6 
Satellites  for  best  GDOP:  9121619 
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Figure  4.14  Skyplot,  Scenario  9 


Table  4.9  Monte  Carlo  Simulation  Results,  Scenario  9 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

26.29 

GDOP 

3.60 

Condition  H 

987.14 

7.41 

7.23 

Error  x 

0.72 

0.74 

0.75 

2.06 

Error  y 

80.04 

5.98 

6.05 

5.97 

6.78 

Error  z 

-67.76 

-3.68 

-3.51 

-0.86 

Error  bias 

344.95 

-3.58 

-4.02 

-2.80 

Miss  distance 

106.81 

7.06 

7.15 

6.97 

7.13 

Estimated  std  x 

770.86 

Estimated  std  y 

2920.93 

Estimated  std  z 

2526.94 

Estimated  std  bias 

13741.87 

Experimental  std  x 

2851.21 

87.03 

87.23 

86.21 

119.07 

Experimental  std  y 

10778.94 

186.99 

186.53 

212.71 

Experimental  std  z 

9322.34 

132.66 

132.24 

208.80 

Experimental  std  bias 

50673.38 

144.76 

152.38 

HI 

Expected  std  x 

2932.02 

80.72 

BH 

Expected  std  y 

11109.87 

207.21 

186.80 

212.20 

Expected  std  z 

9611.34 

206.77 

Expected  std  bias 

52267.75 

87.99 

141.49 

165.36 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time;  t  =  81900  sec 
User  Position;  40  N,  105  W,  300  m 
Satellite  Availability;  N  =  6 
Satellites  for  best  GDOP;  15 18 19  21 
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Figure  4.15 


Skyplot,  Scenario  10 


Table  4.10  Monte  Carlo  Simulation  Results,  Scenario  10 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

26.09 

GDOP 

3.26 

3.92 

Condition  H 

490.73 

9.02 

Error  x 

-6.54 

0.05 

0.02 

0.17 

Error  y 

-25.03 

0.06 

0.09 

-0.34 

Error  z 

22.86 

0.88 

1.06 

-0.14 

Error  bias 

-115.74 

0.79 

0.63 

0.14 

Miss  distance 

34.52 

0.89 

0.88 

1.06 

Estimated  std  x 

538.12 

Estimated  std  y 

2048.11 

Estimated  std  z 

1795.61 

Estimated  std  bias 

9513.75 

Experimental  std  x 

2027.26 

77.76 

77.88 

77.36 

89.65 

Experimental  std  y 

7697.27 

147.92 

146.93 

205.44 

Experimental  std  z 

6736.57 

215.51 

215.53 

215.25 

250.14 

Experimental  std  bias 

35719.62 

177.29 

180.92 

175.89 

197.01 

Expected  std  x 

2062.54 

68.85 

77.84 

89.96 

Expected  std  y 

7850.15 

■Ebkm 

207.83 

Expected  std  z 

6882.36 

251.05 

Expected  std  bias 

36465.09 

124.86 

198.61 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time:  t  =  74700  sec 
User  Position:  40  N,  105  W,  300  m 
Satellite  Availability:  N  =  6 
Satellites  for  best  GDOP:  11 12  16 18 


Figure  4.16  Skyplot,  Scenario  11 


Table  4.11  Monte  Carlo  Simulation  Results,  Scenario  11 


Closed 

Form 

Maximum 

Likelihood 

Kalman 

Update 

Iterative 

Least 

Squares 

Best 

Four 

Estimated  a 

39.08 

GDOP 

3.97 

4.39 

Condition  H 

2442.63 

11.34 

9.93 

Error  x 

-145.85 

-0.72 

-0.72 

-0.92 

Error  y 

-540.48 

-5.62 

-5.68 

-5.34 

-7.86 

Error  z 

467.22 

5.96 

5.99 

5.70 

maimm 

Error  bias 

-2500.71 

4.84 

4.51 

5.52 

Miss  distance 

729.17 

8.22 

8.29 

7.84 

11.42 

Estimated  std  x 

2540.03 

Estimated  std  y 

9358.44 

Estimated  std  z 

8071.08 

Estimated  std  bias 

43840.47 

Experimental  std  x 

6472.81 

97.25 

97.26 

97.25 

Experimental  std  y 

23864.01 

175.63 

176.13 

Experimental  std  z 

20581.58 

263.56 

Experimental  std  bias 

111746.08 

223.90 

227.19 

218.73 

Expected  std  x 

6499.87 

97.35 

96.90 

108.92 

Expected  std  y 

23948.04 

174.18 

206.69 

Expected  std  z 

20653.72 

265.02 

290.10 

Expected  std  bias 

112186.77 

178.45 

218.69 

231.89 
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EXPERIMENTAL  SCENARIO  DATA 

GPS  time;  t  =  89100  sec 
User  Position;  40  N,  105  W,  300  m 
Satellite  Availability;  N  =  5 
Satellites  for  best  GDOP;  12 15 17  21 
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Figure  4.17 


Skyplot,  Scenario  12 


Table  4.12  Monte  Carlo  Simulation  Results,  Scenaxio  12 


Iterative 

Closed 

Maximum 

Kalman 

Least 

Best 

Form 

Likelihood 

Update 

Squares 

Four 

Estimated  a 

GDOP 

3.89 

Condition  H 

220.92 

9.81 

10.54 

Error  x 

-9.38 

-2.49 

-2.48 

-2.48 

Error  y 

-26.46 

Error  z 

22.53 

0.05 

Error  bias 

-123.96 

-0.89 

-1.05 

-0.73 

-1.21 

Miss  distance 

36.00 

2.49 

2.50 

Estimated  std  x 

Estimated  std  y 

Estimated  std  z 

Estimated  std  bias 

Experimental  std  x 

1333.94 

102.45 

102.57 

102.43 

102.43 

Experimental  std  y 

5128.78 

290.92 

292.34 

289.14 

366.57 

Experimental  std  z 

4356.99 

134.99 

132.57 

138.12 

Experimental  std  bias 

23750.81 

204.84 

221.34 

201.19 

230.79 

Expected  std  x 

1344.70 

103.09 

101.49 

101.51 

Expected  std  y 

5146.84 

Expected  std  z 

4374.34 

126.54 

131.95 

137.47 

Expected  std  bias 

23947.62 

281.68 

200.68 

230.91 
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4.11  Experimental  Results,  Ground-Based  Pseudolite  Test  Environment 

The  Monte  Carlo  simulation  results  associated  with  the  ground-based  pseudo¬ 
lite  test  environment  are  presented  in  Table  4.13.  This  test  environment  includes 
three  test  scenarios  A,  B,  and  C  as  described  in  Section  4.3. 

4.11.1  Analysis  of  Results,  Conventional  Iterative  Least  Squares  Algorithm. 

The  results  for  the  conventional  ILS  algorithm  are  discussed  first,  as  they  will  serve 
as  a  baseline  for  comparison.  The  results  in  Table  4.13  show  that  the  ILS  algorithm 
is  capable  of  producing  fairly  good  estimates  of  the  four  GPS  estimation  parameters. 
The  results  also  show  that  the  accuracy  of  the  z  user  position  coordinate  estimate 
degrades  as  the  GPS-like  transmitter  or  pseudolite  below  the  user  position  is  moved 
away  from  directly  below  the  user.  In  this  given  experimental  test  environment,  with 
the  user  at  an  altitude  of  10000  m,  it  was  observed  experimentally  that  there  must 
be  a  pseudolite  within  approximately  500  m  from  the  user’s  projection  on  ground 
in  order  to  produce  a  geometry  that  can  be  handled  by  the  iterative  algorithm.  As 
the  user  altitude  decreases,  the  problem  is  amplified  and  a  smaller  offset  can  be 
tolerated.  This  problem  associated  with  using  the  ILS  algorithm  in  ground-based 
planar  arrays  is  known  and  led  to  the  following  recommendations  by  McKay  [20] : 

•  The  radius  of  the  outer  pattern  must  be  as  large  as  possible,  essentially  at  the 
horizon. 

•  The  user  altitude  must  be  kept  as  large  as  possible  essentially  at  its  flight 
envelope  limit. 

•  The  user  must  be  flown  along  a  known  ground  array  test  pattern  which  guar¬ 
antees  that  a  pseudolite  is  always  close  to  being  directly  below  the  user. 
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Table  4.13  Monte  Carlo  Simulation  Results,  Ground  Test  Environment 


Closed  Form 


Scenario 


Estimated  a 
condition  H 


Error  x 
Error  y 


Error  z 


Error  bias 
Miss  distance 


Estimated  std  x 
Estimated  std  y 


Estimated  std  z 


Estimated  std  bias 
Experimental  std  x 


Experimental  std  y 
Experimental  std  z 


Experimental  std  bias 


Kalman  Update  |  Error  x _ 

Error  y 


Error  z 
Error  bias 


Miss  distance 


Estimated  std  x 
Estimated  std  y 


Estimated  std  z 
Estimated  std  bias 


Experimental  std  x 


Experimental  std  y 
Experimental  std  z 


Experimental  std  bias 
Error  x 


Error  y 


Error  z 
Error  bias 


Miss  distance 
Experimental  std  x 


Experimental  std  y 


Experimental  std  z 
Experimental  std  bias 


A 


99.86 

953483.12 


B 

C 

95.56 

95.22 

1216172.19 

1204038.42 

-0.53 

0.12 

0.64 

-0.14 

293416.35 

1113025.86 

264.20 

259.60 

293416.35 

1113025.86 

28.81 

28.73 

28.81 

28.72 

25016822.93 

24938973.83 

295.77 

296.26 

33.17 

33.15 

1309728.17 


260.13 

1309728.17 


30.65 

30.17 


19937827.70 


312.19 

33.18 


251.55 


9.85 

56.98 


-9632.47 

312.34 


9632.65 


542.53 


3086.50 


mm 
EKyihPia 


36.80 


124.68 

498.92 


368.13 


235.98 


-0.52 

0.44 


-10151.91 

264.59 


10151.91 


33.27 


37.84 


9904.77 

238.59 


28.89 


33.05 

367.70 


239.60 


-9872.20 

260.17 


9872.20 


0.97 

33.73 


32.61 


341.74 

244.70 


9903.54 

241.33 


28.79 


0.27 


0.19 


11.46 

-10.62 


11.47 

33.71 


33.91 


347.54 

249.12 
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The  results  for  Scenario  C  show  the  case  when  the  pseudolite  below  the  user  is 
offset  by  400  m,  which  is  close  to  the  largest  offset  that  can  be  tolerated  without  pro¬ 
ducing  a  singular  geometry.  In  a  singular  geometry,  a  solution  cannot  be  obtained, 
as  was  the  case  in  Scenario  A  when  no  pseudolite  was  available  directly  below  the 
receiver  or  in  the  vicinity.  Scenario  B  whose  layout  is  shown  in  Figure  4.1,  in  Se- 
cion  4.3,  represents  an  optimum  geometry;  consequently,  the  estimated  parameter 
errors  were  all  small.  The  results  for  Scenario  C  show  little  change  on  the  estimates 
of  the  X  and  y  user  position  coordinates  and  the  range-equivalent  user  clock  bias, 
but  the  error  in  the  jj;  user  position  coordinate  estimate  increases  form  0.39  m  to 
11.46  m.  The  standard  deviation  on  the  z  user  position  coordinate  estimation  error 
was  one  order  of  magnitude  larger  than  for  the  x  and  y  coordinate;  however,  there 
was  no  impact  on  the  2;  user  position  coordinate  error  standard  deviation  caused  by 
the  400  m  offset  in  the  pseudolite  below  the  user. 

4.11.2  Analysis  of  Results,  Closed-Form  Algorithm.  The  results  in  Ta¬ 
ble  4.13  indicate  that,  in  this  test  environment,  the  closed-form  algorithm  produces 
excellent  estimates  of  the  pseudorange  measurement  noise  strength  cr.  The  a  pro¬ 
duced  by  the  algorithm  ranged  from  95.22  m  to  99.86  m  compared  with  the  true 
value  of  100  m  used  in  simulating  the  measurement  noise.  This  test  environment 
was  intentionally  simulated  using  ground-based  planar  arrays  of  36  or  37  pseudolites 
because  from  the  experimental  results  presented  in  Figure  4.5,  in  Section  4.8.2,  it 
was  anticipated  that  this  number  of  transmitters  would  be  required  to  obtain  reli¬ 
able  estimates  of  a.  The  pseudolite  availability  in  this  test  environment  is  greater 
than  the  satellite  availability  that  can  be  achieved  in  conventional  NAVSTAR  GPS 
scenarios.  High  availability  is  favorable  for  the  estimation  of  a  since  the  validity  of 
the  ergodic  assumption  invoked  in  Equation  (3.20),  which  formed  the  basis  of  the 
estimation  of  cr,  is  strengthened  by  increasing  transmitter  availability.  Furthermore, 
Scenario  A  yielded  the  most  favorable  conditions  for  estimating  the  measurement 
noise  strength  since  the  pseudoranges  to  all  transmitters  were  essentially  the  same. 
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validating  the  average  pseudorange  approximation  used  in  Equation  (3.5).  Unlike 
the  results  obtained  for  the  simulated  GPS  test  environment,  the  estimation  of  Ux 
and  Uy  user  position  coordinates  are  extremely  good,  with  errors  marginally  smaller 
than  those  obtained  with  the  conventional  ILS  algorithm;  however,  the  user  posi¬ 
tion  coordinate  estimation  error  is  outrageously  large,  ranging  from  2.9  x  10®  meters 
to  1.3  X  10®  meters.  It  is  evident  from  the  results  that  the  closed-form  algorithm 
has  a  serious  observability  problem  along  the  2;  coordinate  axis.  The  problem  is 
marginally  worst  in  Scenario  A  where  there  is  no  transmitter  below  the  user  along 
the  2r  coordinate  axis.  The  miss  distance  results  are  not  used  for  comparison  pur¬ 
pose  because  they  reflect  the  overly  dominant  error  in  the  z  coordinate  estimates. 
Given  the  algorithm’s  ability  to  estimate  a  accurately,  the  estimation  of  the  stan¬ 
dard  deviation  of  the  estimation  error  is  also  extremely  good,  where  the  estimates 
were  all  within  20.24  percent  of  the  experimentally  derived  standard  deviations  for 
all  four  estimation  parameters,  in  the  three  scenarios.  For  the  x  and  y  coordinates, 
the  experimental  standard  deviation  is  slightly  lower  than  that  achieved  with  the 
ILS  algorithm.  The  standard  deviation  on  the  bias  estimation  error  is  slightly  larger 
than  that  achieved  with  the  ILS  algorithm.  Whether  or  not  a  transmitter  is  di¬ 
rectly  below  the  user,  had  virtually  no  effect  on  the  effect  on  the  performance  of 
the  closed-form  algorithm.  The  closed-form  algorithm  is  not  dependent  on  having  a 
transmitter  directly  below  at  all  times,  which  would  be  physically  impossible.  The 
closed-form  algorithm  is  capable  of  producing  x  and  y  position  coordinate  estimates 
better  than  the  conventional  iterative  algorithm;  however,  it  does  significantly  worse 
at  estimating  the  2:  position  coordinate  and  slightly  worse  at  estimating  the  range- 
equivalent  user  clock  bias.  In  addition,  given  the  extremely  low  estimation  errors  in 
the  X  and  y  user  position  coordinate  estimates,  it  appears  that  the  geometry  pro¬ 
duced  by  a  spread-out  ground-based  planar  array  of  pseudolites  is  more  favorable 
to  the  closed-form  algorithm  than  any  geometry  that  can  be  produced  considering 
strictly  the  24  satellite  NAVSTAR  GPS  constellation.  If  it  was  necessary  to  estimate 
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only  the  x  and  y  user  position  coordinates  from  signals  obtained  from  pseudolites  in 
a  ground-based  planar  array,  this  c/osed-/orm  . algorithm  would  be  the  algorithm  of 
choice. 

4.11.3  Analysis  of  Results,  Kalman  Update.  The  results  indicate  that, 
starting  from  the  four  GPS  parameter  estimates  produced  by  the  closed-form  al¬ 
gorithm  and  applying  the  Kalman  update  algorithm  twice,  there  is  a  potential  for 
marginal  improvement  or  corruption  of  the  estimates,  of  the  x  and  y  user  position 
coordinates  and  the  range-equivalent  user  clock  bias.  The  small  changes  in  the  esti¬ 
mates  of  the  X  and  y  user  position  coordinates  allowed  the  Kalman  update  algorithm 
to  produce  greatly  improved  estimate  of  the  2  user  position  coordinate.  The  stan¬ 
dard  deviation  of  x  and  y  user  position  coordinates  and  the  range-equivalent  user 
clock  bias  estimates  are  also  effected  very  little,  sometimes  for  the  worse,  other  times 
for  the  better.  The  algorithm’s  capability  to  estimate  the  standard  deviation  on  the 
X  and  y  user  position  coordinates  and  the  range  equivalent  user  clock  bias  also  works 
well.  The  estimation  errors  for  the  2  user  position  coordinate  range  from  9632.47  m 
to  10151.91  m,  representing  two  orders  of  magnitude  in  error  reduction  compared  to 
the  closed-form  algorithm.  Despite  this  improvement,  the  error  in  the  2  coordinate  is 
still  too  large  to  render  this  algorithm  useful.  The  improvement  in  the  experimental 
standard  deviation  in  the  2  coordinate  estimate  is  more  significant  than  the  improved 
estimate  itself;  however,  the  algorithm  is  grossly  overestimating  its  error  on  the  2 
coordinate  estimate.  Based  on  the  results,  the  Kalman  update  algorithm  does  not 
prove  very  useful  in  this  ground-based  planar  array  test  environment,  as  it  does  not 
provide  any  significant  improvement  over  the  closed-form  algorithm;  moreover,  the 
risk  of  corrupting  the  x  and  y  user  position  coordinate  estimates  exists. 

4.12  Summary  of  Experimental  Results 

In  typical  near-earth  GPS  scenarios,  the  root  mean  square  miss  distances  as¬ 
sociated  with  the  closed-form  algorithm,  ranging  from  4.93  m  to  729.17  m,  were 
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much  worse  than  the  baseline  ILS  algorithm’s  results,  which  produced  a  maximum 
miss  distance  of  7.84  m.  However,  when  the  closed-form  algorithm  is  assisted  by 
a  supplementary  algorithm,  the  miss  distances  were  reduced  to  the  baseline  level, 
within  8.22  m  and  8.29  m  for  the  maximum  likelihood  and  Kalman  update  algorithms 
respectively.  Based  only  on  miss  distance,  the  deficiency  in  the  unassisted  closed- 
form  algorithm  is  evident;  however,  the  experimental  standard  deviations  associated 
with  the  four  GPS  solution  parameters,  which  are  much  more  reliable  performance 
indicators,  clearly  reveal  the  seriousness  of  the  problem.  The  experimental  standard 
deviation,  on  the  three  position  parameter  estimates,  yielded  by  the  closed-form 
algorithm  ranged  from  1078.08  m  to  23864.01  m,  compared  to  the  baseline  ILS 
results,  which  ranged  from  55.56  m  to  289.14  m.  The  experimental  standard  devi¬ 
ation  on  the  range-equivalent  user  clock  bias  estimates  were  much  worse,  ranging 
from  19014.42  m  to  111746.08  m.  The  maximum  likelihood  and  Kalman  update  al¬ 
gorithms  yielded  experimental  standard  deviations  comparable  to  the  baseline  ILS 
results  for  the  three  position  parameter  estimates  but  consistently  slightly^larger 
than  the  baseline  ILS  results  for  the  user  clock  bias.  The  closed-form  algorithm’s 
performance  is  inadequate  for  stand-alone  use;  however,  its  results  are  used  as  ini¬ 
tialization  for  the  supplementary  algorithms,  to  achieve  performance  comparable  to 
that  of  the  baseline  ILS  algorithm.  The  estimation  error  covariance  prediction  per- 
formance  of  the  closed-form  algorithm  is  dependent  on  how  well  the  pseudorange 
measurement  noise  variance  (cr^)  can  be  estimated,  which  is  a  function  of  satellite 
availability.  Under  typical  satellite  availability  conditions  it  was  not  possible  to  re¬ 
trieve  good  data-driven  estimates  of  cr,  but  as  satellite  availability  is  increased,  the 
estimates  of  a  improve.  Using  a  scenario  concatenation  scheme,  a  mean  estimated 
a  of  92.5  m  was  retrieved,  which  compares  well  with  the  true  value  of  100  m  used  to 
generate  the  experimental  measurement  noise. 

The  unconventional  geometries  in  the  ground-based  planar  array  scenarios  us¬ 
ing  pseudolites,  produced  results  that  were  more  favorable  to  the  closed-form  al- 


4-43 


gorithm.  The  horizontal  positioning  performance  of  the  closed-form  algorithm  was 
marginally  better  than  the  baseline  ILS  results  both  in  terms  of  mean  errors  and 
experimental  standard  deviation.  The  biggest  attribute  of  the  closed-form  algorithm 
was  that,  unlike  the  ILS  algorithm  which  is  highly  dependent  on  having  a  pseudo- 
lite  directly  below  the  user,  the  closed-form  algorithm  can  produce  good  horizontal 
positioning  results  under  any  geometry.  Given  the  large  pseudolite  availability  in 
the  ground-based  planar  array,  excellent  estimates  a  cr,  between  95.2  m  and  99.9  m, 
were  retrieved.  As  a  result,  the  closed- form  algorithm’s  error  covariance  prediction 
performance  was  very  good. 


4-44 


V.  Conclusions  and  Recommendations 

This  chapter  presents  a  brief  summary  of  the  performance  related  issues  of 
the  closed-form  algorithm  and  the  supporting  algorithms  developed  for  this  thesis 
research.  Emphasis  is  placed  on  the  identifying  the  areas  of  strength  for  the  algo¬ 
rithms  and  suggesting  applications  for  which  they  are  best  suited.  The  shortcomings 
seen  in  the  Chapter  IV  related  to  algorithm  deficiency  and  inconclusive  findings  are 
discussed,  and  potential  approaches  to  correcting  the  problems  are  presented.  The 
chapter  sums  up  with  recommendations  for  future  work  that  can  prove  rewarding  in 
bringing  the  concept  of  using  a  stochastic  closed-form  GPS  position  determination 
algorithm  closer  to  reality. 

5.1  Typical  Near-Earth  GPS  Applications 

When  compared  to  the  current  linear  iterative  algorithms  using  approximative 
linearization,  the  use  of  stochastic  modeling  and  estimation  techniques  applied  to 
GPS,  describes  a  more  correct  approach  to  solving  the  system  of  nonlinear  GPS 
pseudorange  equations  from  a  mathematical  perspective.  Although  the  concept 
is  sound,  the  results  yielded  by  this  research  indicate  that  there  are  no  benefits 
to  be  gained  in  terms  of  estimation  accuracy  when  dealing  with  the  most  typical 
GPS  application  of  near-earth  navigation  using  strictly  the  satellites  available  in  the 
NAVSTAR  GPS  24  satellite  constellation.  Under  near-earth,  large  SNR  conditions, 
the  research  demonstrates  that  the  stand-alone  closed-form  algorithm  is  incapable  of 
producing  a  reasonable  GPS  solution  estimate  without  being  supplemented  by  either 
the  maximum  likelihood  or  the  Kalman  update  algorithm.  The  performance  of  both 
supplementary  algorithms  is  essentially  equivalent  but  the  latter  is  the  preferred  al¬ 
gorithm  because  it  is  closed-form  in  nature  and  because  it  possesses  the  capability 
of  predicting  its  estimation  error  covariance.  For  these  reasons,  the  conclusions  will 
emphasize  results  associated  with  the  Kalman  update  algorithm. 
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This  research  clearly  indicates  that  the  developed  algorithm,  consisting  of  the 
closed-form  algorithm  as  supplemented  by  the  Kalman  update  algorithm,  does  not 
offer  any  significant  performance  improvements  over  the  conventional  ILS  algorithm 
in  typical  near-earth  navigation  applications.  The  relative  performance  of  the  devel¬ 
oped  algorithm  can  be  summed  up  by  the  following  points: 

1.  The  performance  in  terms  of  accuracy,  based  on  mean  miss  distance,  of  the 
developed  algorithm  is  equivalent  to  that  of  the  conventional  ILS  algorithm. 
In  all  the  realistic  near-earth  test  scenarios  analyzed  in  this  thesis,  the  miss 
distances  yielded  by  the  developed  algorithm  were  a  fraction  of  a  meter  either 
above  or  below  those  yielded  by  the  ILS  algorithm. 

2.  The  performance  in  terms  of  experimental  standard  deviation  of  the  GPS  pa¬ 
rameter  estimates  of  the  developed  algorithm  is,  at  best,  slightly  degraded  when 
compared  to  that  of  the  conventional  ILS  algorithm.  This  performance  indica¬ 
tor  is  a  measure  of  how  sensitive  the  algorithm  is  to  pseudorange  measurement 
noise  and  establishes  the  estimation  error  bounds.  The  experimental  standard 
deviation  associated  with  the  three  position  parameter  estimates  {ux^Uy^uf) 
produced  by  the  developed  algorithm  where  marginally  larger;  all  within  4.6 
percent  above  the  ILS  results.  The  performance  degradation  of  the  experimen¬ 
tal  standard  deviation  associated  with  the  estimate  of  the  range-equivalent 
user  clock  bias  was  more  significant,  ranging  from  2.8  to  35.3  percent  above 
the  ILS  results. 

3.  Considering  only  satellites  available  in  the  NAVSTAR  GPS  constellation,  there 
is  insufficient  satellite  availability,  preventing  the  closed-form  algorithm  from 
producing  a  reliable  estimate  of  the  pseudorange  measurement  noise  strength. 
The  research  showed  improvement  in  estimating  the  noise  strength,  cr,  with 
increasing  satellite  availability;  however,  an  availability  of  at  least  36  satellites 
is  needed  to  produce  a  reliable  estimate.  Since  this  number  is  not  achievable  in 
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typical  near-earth  GPS  applications,  the  developed  algorithm  cannot  predict 
a  reliable  estimation  error  covariance,  under  these  conditions. 

4.  When  the  algorithms  are  artificially  informed  of  the  pseudorange  measurement 
noise  strength,  cr,  the  developed  algorithm  as  well  as  the  conventional  ILS  algo¬ 
rithm  performed  extremely  well  at  producing  the  estimation  error  covariance 
matrix.  The  diagonal  elements  of  the  covariance  matrix  were  used  to  calculate 
the  expected  standard  deviation  of  the  individual  GPS  parameter  estimates 
for  comparison  to  the  experimentally  derived  values.  The  developed  algorithm 
performed  slightly  worse  than  the  conventional  ILS  algorithm.  The  expected 
standard  deviations  for  the  position  parameters  for  the  developed  algorithm 
were  all  within  14.0  percent  of  the  experimental  values  and  the  standard  devi¬ 
ations  on  the  clock  bias  were  off  by  as  much  as  44.2  percent.  These  results  do 
not  compare  very  well  to  the  maximum  of  2.25  percent  on  all  four  parameters 
achieved  by  the  ILS  algorithm. 

The  largest  deficiency  for  the  developed  algorithm  is  associated  with  the  estimation 
of  the  clock  bias,  which  is  not  an  issue  if  only  positioning  accuracy  is  considered.  It 
is  suspected  that  this  deficiency  is  a  direct  result  of  the  composition  of  the  regressor 
matrix.  The  closed-form  regressor  matrix  is  such  that  the  last  column  is  made  up 
of  differenced  pseudorange  measurements,  whereas  the  the  first  three  columns  are 
made  up  of  differenced  satellite  position  coordinates.  This  composition  scheme  leads 
to  an  ill-conditioned  regressor  that  is  unfavorable  to  the  estimation  of  the  last  GPS 
parameter,  namely  the  range-equivalent  user  clock  bias. 

5.2  Unconventional  GPS  Applications 

The  developed  algorithm  which  consists  of  the  basic  closed-form  algorithm  as 
supplemented  by  the  Kalman  update  algorithm  is,  as  a  whole,  a  closed-form  algo¬ 
rithm.  Inherently,  it  does  not  require  any  form  of  user  initialization  in  order  to 
produce  a  GPS  solution.  One  direct  benefit  of  a  closed-form  GPS  position  deter- 
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mination  algorithm  is  the  ability  to  produce  a  solution  in  geometrical  conditions 
very  different  from  the  typical  near-earth  applications  for  which  GPS  is  known  to 
be  best  suited.  These  uncommon  geometrical  conditions  can  be  encountered  in  deep 
space  applications  and  in  pseudolite  ground-based  arrays,  both  of  which  can  lead 
to  difficulties  for  the  conventional  iterative  algorithms.  In  these  applications,  the 
developed  algorithm  alleviates  any  concern  of  multiple  solutions  or  of  convergence 
towards  the  wrong  solution  as  a  result  of  incorrect  initialization,  as  is  the  case  for 
the  conventional  iterative  algorithms. 

This  thesis  research  did  not  include  any  deep  space  simulated  scenarios  for 
testing  purpose;  hence,  no  comments  can  be  made  on  performance  issues  in  deep 
space.  However,  the  developed  algorithm  was  tested  in  simulated  pseudolite  ground- 
based  array  scenarios  and  it  was  found  to  be  much  more  tolerant  to  geometrical 
conditions  than  was  the  conventional  ILS  algorithm.  The  developed  algorithm  is  not 
dependent  on  having  a  pseudolite  in  the  vicinity  of  user’s  position  projection  on  the 
ground.  The  conventional  ILS  algorithm  was  incapable  of  producing  a  solution  if  a 
pseudolite  was  not  present  virtually  directly  below  the  user  position.  The  acceptable 
offset  tolerance  for  the  receiver  directly  below  the  user  is  proportional  to  the  ground 
array  size  and  the  user  altitude.  In  a  small  size  ground  array  with  the  user  at  low 
altitudes,  the  tolerance  would  be  extremely  small,  requiring  that  a  pseudolite  always 
be  precisely  below  the  user.  This  condition  is  extremely  restrictive  on  the  flight 
pattern  and  would  require  pseudolite  placed  in  close  proximity  to  each  others  below 
the  intended  flight  pattern.  All  these  concerns  are  alleviated  by  using  the  developed 
algorithm. 

In  estimating  strictly  the  horizontal  user  position  coordinates  Ux  and  Uy,  the 
closed-form  algorithm  on  its  own  without  the  assistance  of  a  supplementary  algo¬ 
rithm  was  the  best  performing  algorithm.  The  mean  error  on  these  two  estimation 
parameters  was  smaller  as  was  the  experimental  standard  standard  deviation,  when 
compared  to  the  results  yielded  by  the  baseline  ILS  algorithm.  The  algorithm’s  in- 
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sensitivity  to  ground  array  geometry  and  its  extremely  good  horizontal  positioning 
accuracy  makes  it  well  suited  for  test  range  applications.  Nonetheless,  the  closed- 
form  algorithm  suffers  from  serious  observability  problems  along  the  vertical  axis, 
as  was  observed  from  the  extremely  poor  accuracy  in  the  estimation  of  vertical  user 
position  coordinate,  The  observability  problem  along  the  vertical  axis  is  a  re¬ 
sult  of  how  the  planar  array  is  constructed.  If  the  array  were  truly  planar,  the  z 
coordinates  of  all  the  pseudolites  would  be  exactly  the  same,  which  would  result 
in  a  singular  regressor  matrix.  Recalling  that  the  third  column  of  the  regressor  is 
made  up  of  differenced  satellite  z  position  coordinates,  in  a  truly  planar  array  the 
third  column  would  be  a  column  of  zeros.  Fortunately,  the  world  is  not  perfectly 
flat  and  it  was  necessary  to  introduce  some  small  random  variations  in  the  simulated 
pseudolite  2;  position  coordinates  to  reflect  terrain  perturbations.  As  the  variations 
in  the  pseudolite  altitudes  is  increased,  the  accuracy  of  the  estimated  vertical  user 
position  coordinate  improves  at  the  expense  of  slightly  degrading  horizontal  posi¬ 
tioning  accuracy.  It  is  probably  best  to  maintain  the  horizontal  positioning  accuracy 
at  the  highest  achievable  level  and  solving  the  vertical  positioning  accuracy  prob¬ 
lem  by  using  an  extremely  accurate  altitude  measuring  instrument  such  as  a  radar 
altimeter  in  conjunction  with  the  closed-form  algorithm  in  the  GPS  receiver.  This 
combination  is  generally  beneficial,  since  GPS,  in  general,  produces  vertical  errors 
that  are  larger  than  the  horizontal  errors. 

5.S  Estimation  Error  Covariance 

One  of  the  driving  motivations  to  this  thesis  research,  in  applying  stochastic 
modeling  and  estimation  techniques  to  the  GPS  pseudorange  equations,  was  the  pos¬ 
sibility  of  producing  an  estimate  of  the  estimation  error  covariance.  This  capability 
would  provide  a  confidence  factor  associated  with  the  GPS  solution  in  real  time.  The 
research  demonstrated  that  it  is  possible  to  estimate  the  pseudorange  measurement 
noise  strength  (<t)  accurately,  using  the  developed  algorithm,  if  satellite  availability  is 
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greater  than  36.  Once  the  value  of  a  is  known,  the  estimation  error  covariance  can  be 
calculated  for  both  the  closed  form  algorithm  developed  in  this  thesis  or  for  the  con¬ 
ventional  iterative  algorithms.  Due  to  the  deterministic  nature  of  the  conventional 
iterative  algorithms,  they  have  no  way  of  estimating  the  pseudorange  measurement 
noise  strength.  Two  ways  of  estimating  cr  were  shown;  one  using  the  concatenation  of 
realistically  simulated  satellite  scenarios  and  the  other  using  simulated  ground-based 
pseudolite  planar  array  scenarios. 

The  concatenation  of  realistic  satellite  scenarios  at  one  hour  intervals  to  obtain 
an  effective  satellite  availability  of  42  proved  to  be  effective  for  the  estimation  of  cr. 
Using  this  approach  while  maintaining  the  user  position  stationary  produced  a  mean 
estimated  cr  of  92.5  meters  which  compares  very  well  to  the  cr  of  100  meters  that  was 
used  to  simulate  the  experimental  noise.  This  approach  is  extremely  restrictive  due 
to  the  fact  that  the  user  position  is  maintained  stationary  at  the  measurement  point, 
therefore,  its  use  for  navigation  in  flight  appears  limited.  Although  the  algorithm 
can’t  be  used  for  estimation  purposes  during  flight  in  real  time,  it  could  be  used  to 
calibrate  the  receiver  if  it  is  known  that  it  will  be  used  in  the  vicinity  of  the  measure¬ 
ment  point.  Receiver  calibration  can  assist  in  improving  the  receiver’s  accuracy  by 
compensating  for  noise  associated  with  the  measurements  as  well  as  the  noise  that 
it  is  responsible  for  in  itself.  Obviously,  the  closer  the  user  is  to  the  measurement 
point,  the  greater  the  benefits  of  the  receiver  calibration. 

The  approach  of  estimating  the  noise  strength  from  a  fixed  measurement  point 
is  very  well  suited  for  GPS  applications  in  which  the  user  position  is  stationary, 
such  as  in  site  surveying.  Site  surveying,  is  becoming  an  increasingly  popular  ap¬ 
plication  for  GPS.  In  site  surveying,  accuracy  is  of  extreme  importance.  The  use 
of  the  measurement  noise  estimation  scheme  will  allow  for  the  calculation  of  the 
estimation  error  covariance  which  can  be  used  as  a  confidence  indicator  as  well  as  a 
means  to  demonstrate  that  a  desired  level  of  accuracy  has  been  achieved.  Once  the 
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measurement  noise  strength  has  been  estimated,  it  can  be  used  within  the  algorithm 
of  choice  to  produce  the  estimation  error  covariance. 

Using  a  ground  array  with  36  pseudolites  in  a  circular  pattern  produced  a 
mean  estimated  a  of  99.86  meters  which  is  remarkably  close  to  the  actual  experi¬ 
mental  noise  strength  of  100  meters.  This  concept  could  prove  extremely  useful  in 
test  range  arrangements  such  as  the  SARS  [25],  [20].  It  has  the  potential  of  being 
useful  in  a  Wide  Area  Augmentation  Systems  (WAAS)  [15]  which  is  due  to  become 
operational  in  the  near  future.  The  WAAS  is  a  wide-area  DGPS  concept  being  de¬ 
veloped  by  the  Federal  Aviation  Administration  (FAA)  to  increase  GPS  availability 
in  support  of  flight  operations  from  en-route  navigation  through  Category  I  preci¬ 
sion  approaches.  WAAS  will  consist  of  a  ground-based  communications  network  and 
several  geosynchronous  satellites  to  provide  global  coverage.  Differential  corrections 
and  additional  ranging  data  will  be  broadcast  to  users  from  geostationary  satellites 
using  a  signal  similar  to  the  NAVSTAR  GPS  signal  [15].  The  integrated  use  of  GPS, 
GLONASS,  and  the  proposed  WAAS  could  form  a  Global  Navigation  Satellite  Sys¬ 
tem  (GNSS)  that  would  guarantee  much  improved  satellite  availability  that  would 
allow  reliable  estimations  of  the  pseudorange  measurement  noise  strength  using  the 
concepts  presented  in  this  thesis.  In  addition  to  the  WAAS,  the  possibility  of  other 
GPS  augmentation  systems  may  be  implemented  for  other  specific  purposes,  further 
improving  satellite  availability. 

5.4  Condition  Number  of  the  Regressor 

For  the  conventional  iterative  algorithms,  this  research  demonstrated  that  there 
exists  a  direct  relationship  between  the  estimation  accuracy  and  the  GDOP.  This 
relationship  extended  very  well  to  the  condition  number  of  the  regressor,  H,  where 
both  measures  were  essentially  equivalent  indicators  of  geometrical  strength.  In  the 
case  of  the  closed-form  algorithms,  although  GDOP  is  defined  mathematically,  it 
has  no  physical  meaning,  since  it  is  formulated  for  use  specifically  with  the  matrix 
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of  direction  cosines  associated  with  iterative  algorithms;  hence,  the  only  suitable 
measure  is  the  condition  number  of  the  regressor.  The  sensitivity  of  the  closed-form 
algorithms  to  geometrical  conditions  is  not  well  known  or  understood  but  it  is  ob¬ 
vious  that  what  constitutes  good  geometry  for  the  conventional  iterative  algorithms 
does  not  necessarily  constitute  good  geometry  for  the  closed-form  algorithm.  A 
loose  relationship  between  the  condition  number  of  the  regressor  and  the  achievable 
estimation  accuracy  was  observed  for  the  closed-form  algorithm.  This  relationship 
cannot  be  generalized  into  an  error  bound  concept,  as  is  the  case  for  iterative  algo¬ 
rithms,  because  the  algorithm’s  behavior  changes  drastically  for  different  scenarios  or 
applications.  The  pseudolite  ground-based  array  scenarios  were  extremely  favorable 
to  the  closed-form  algorithms,  in  the  estimation  of  horizontal  user  position  coordi¬ 
nates,  despite  the  extremely  large  condition  number  associated  with  the  regressor 
matrix.  The  manifestation  of  the  ill-conditioned  regressor  was  obvious  from  the  poor 
accuracy  in  estimating  the  range-equivalent  user  clock  bias. 

5.5  Recommendations 

This  section  presents  areas  that  remain  to  be  explored  that  can  be  taken  on  as 
follow  on  research.  The  results  produced  by  this  research,  although  not  very  favorable 
for  the  implementation  of  the  stochastic  closed-form  algorithm,  has  brought  to  light 
many  interesting  issues  that  deserve  further  investigation.  Given  that  this  area  of 
research  is  relatively  new  when  compared  to  the  iterative  algorithms,  further  research 
is  required  to  assess  the  viability  of  applying  stochastic  mathematical  modeling  to 
solve  the  system  of  GPS  pseudorange  equations. 

5.5.1  GPS  Measurement  Noise  Levels  Investigation.  In  view  of  the  out¬ 
come  of  this  research,  the  most  intriguing  dilemma  that  remains  unsolved  is,  why  the 
conventional  ILS  algorithm  performs  so  well  despite  the  approximative  nature  of  the 
linearization  used  and  the  deterministic  approach  taken,  both  of  which  are  sources 
of  error.  This  dilemma  is  more  intuitively  solved  by  rationalizing  why  the  stochastic 
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closed-form  algorithm  in  this  thesis  performed  more  poorly  than  anticipated.  The 
explanation  to  this  dilemma  is  related  to  the  level  of  noise  that  is  being  consid¬ 
ered.  The  low  noise  levels  considered  in  this  research,  which  are  typical  of  GPS, 
are  too  small  to  gain  any  benefit  from  a  weighted  minimum  variance  estimate.  The 
inverse  covariance  weighting  actually  introduced  weighting  errors  into  the  stochastic 
estimate.  Preliminary  testing  indicated  that,  an  unweighted  least  squares  estimate 
may  produce  marginally  better  results  at  the  expense  of  losing  the  capability  of 
producing  the  estimation  error  covariance.  Under  degraded  operating  conditions, 
with  larger  pseudorange  measurement  noise  strengths,  stochastic  estimation  may 
lead  to  better  solutions,  in  terms  of  both  mean  estimation  error  and  experimental 
standard  deviation,  than  the  conventional  iterative  approaches.  Research  is  required 
to  determine  the  noise  levels  at  which  conditions  become  favorable  to  the  stochastic 
modeling  approach  and,  more  importantly,  determining  if  those  noise  levels  can  ever 
be  encountered  with  GPS  in  its  current  state  of  technology. 

5.5.2  Alternate  Stochastic  Closed-Form  Algorithms.  The  derivation  of 
the  two  closed-form  algorithms  in  this  thesis  make  use  of  an  average  pseudorange 
approximation  in  order  to  simplify  the  noise  statistics.  The  approximation  is  also 
needed  to  obtain  an  expression  for  the  covariance  that  is  not  dependent  on  any  of 
the  four  GPS  parameters  being  estimated.  This  approach  enabled  the  algorithm  to 
maintain  a  true  closed-form  structure.  The  derivation  of  the  algorithm  needs  to  be 
revisited  without  the  use  of  any  approximations  in  order  to  evaluate  the  impact  of 
using  the  approximation  through  experimental  analysis.  In  addition,  new  approaches 
to  deriving  alternate  stochastic  closed-form  solutions  to  the  system  of  pseudorange 
equations  must  be  investigated  in  an  attempt  to  obtain  an  algorithm  that  possesses 
the  following  qualities; 

•  The  regressor  matrix  should  have  a  low  condition  number  to  maintain  the 

estimation  error  amplification  bounds  to  a  minimum. 
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•  The  algorithm  should  be  capable  of  producing  an  estimate  of  the  four  GPS 
estimation  parameters  using  only  four  pseudorange  measurements,  as  is  the 
case  for  deterministic  solutions  both  iterative  and  closed-form. 

•  The  algorithm  should  be  capable  of  producing  an  accurate  GPS  solution  with 
a  single  application  without  the  use  of  a  supplementary  algorithm. 

It  must  be  noted  that  the  existence  of,  or  feasibility  of  developing,  an  algorithm  that 
possesses  all  or  any  of  the  above  qualities  is  not  guaranteed. 

5.5.3  Computational  Effectiveness.  Another  area  that  remains  to  be  ex¬ 
plored  is  computational  effectiveness  of  the  closed-form  algorithms.  Given  current 
state  of  computational  power,  this  issue  is  not  one  of  high  priority.  Nonetheless, 
consideration  must  be  given  to  the  computational  effectiveness  to  assess  the  viabil¬ 
ity  of  using  the  algorithms  in  real-time  applications.  The  closed-form  algorithms  are 
expected  to  be  more  effective  from  a  computational  standpoint  than  the  iterative 
algorithms  due  to  the  non-recursive  nature;  however,  the  rewards  of  the  closed-form 
algorithm  may  be  offset  by  an  overly  complex  algorithm  that  requires  computation¬ 
ally  demanding  mathematical  operations. 
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